1!!!============== Find critical point of real space functions 2subroutine topomain 3use topo 4use defvar 5use util 6use function 7implicit real*8(a-h,o-z) 8character c200*200,c1000*1000,ctmp1*20,ctmp2*20,ctmp3*20,ctmp4*20,icp1text*12,icp2text*12 9real*8,allocatable :: randptx(:),randpty(:),randptz(:) !x,y,z of the points in the sphere 10real*8,allocatable :: bassurpathtmp(:,:,:,:) !Used for temporary store bassurpath 11real*8,allocatable :: shanCPrho(:) !For calculate shannon aromaticity 12integer :: shanCPind(100),searchcenlist(500) 13real*8 :: hesstmp(3,3),gradtmp(3) !Temporarily used for calculating curvature 14 15!Initialize searching and plotting parameters 16toposphrad=3D0 !Radius of searching sphere is 3 Bohr 17numsearchpt=1000 !1000 points randomly scattered in the sphere 18sphcenx=0D0 !Position of the sphere center 19sphceny=0D0 20sphcenz=0D0 21!Initialize plot parameter 22ZVU=5.0D0 !Closer than other case 23idrawisosur=0 24ishow3n3=1 25ishow3n1=1 26ishow3p1=1 27ishow3p3=1 28idrawpath=1 29bondradius=0.07D0 !For showing CPs, we use thiner default bond to avoid overlay the (3,-1) 30ratioatmsphere=0.6D0 !Use smaller radius of atom 31textheigh=36 32ishowCPlab=0 33ishowatmlab=0 34ishowpathlab=0 35ishowsearchlevel=0 36ishowattlab=0 !Don't show the result of basin analysis 37 38write(*,*) " !!! Note: All length units in this module are Bohr !!!" 39 40do while(.true.) 41 write(*,*) 42 write(*,"(a)") " ================ Topology analysis ===============" 43 write(*,"(a,i3)") " -11 Delete results and reselect real space function, current:",ifunctopo 44 write(*,*) "-10 Return" 45 write(*,*) "-9 Measure distances, angles and dihedral angles between CPs or atoms" 46 write(*,*) "-5 Modify or print detail or export paths, or calculate property along a path" 47 write(*,*) "-4 Modify or export CPs (critical points)" 48 write(*,*) "-3 Set interbasin surface generating parameters" 49 write(*,*) "-2 Set path searching parameters" 50 write(*,*) "-1 Set CP searching parameters" 51 write(*,*) "0 Print and visualize all generated CPs, paths and surfaces" 52 write(*,*) "1 Search CPs from a given starting point" 53 write(*,*) "2 Search CPs from nuclear positions" 54 write(*,*) "3 Search CPs from midpoint of atom pairs" 55 write(*,*) "4 Search CPs from triangle center of three atoms" 56 write(*,*) "5 Search CPs from pyramid center of four atoms" 57 write(*,*) "6 Search CPs from a batch of points within a sphere" !cube, random in sphere 58 write(*,*) "7 Show real space function values at specific CP or all CPs" 59 write(*,*) "8 Generate the path connected (3,-3) and (3,-1)" 60 write(*,*) "9 Generate the path connected (3,+1) and (3,+3)" 61 write(*,*) "10 Add or delete interbasin surfaces" 62 if (ifunctopo==1) write(*,*) "20 Calculate Shannon aromaticity index" 63 write(*,*) "21 Calculate density curvature perpendicular to a specific plane at a point" 64 read(*,*) isel 65 66 if (isel==-11) then 67 write(*,*) "0 Return" 68 write(*,*) "1 Electron density (Analytical Hessian)" 69 write(*,*) "3 Laplacian of electron density" 70 write(*,*) "4 Value of orbital wavefunction" 71 if (ELFLOL_type==0) write(*,*) "9 Electron localization function (ELF)" 72 if (ELFLOL_type==1) write(*,*) "9 Electron localization function (ELF) defined by Tsirelson" 73 if (ELFLOL_type==2) write(*,*) "9 Electron localization function (ELF) defined by Lu, Tian" 74 if (ELFLOL_type==0) write(*,*) "10 Localized orbital locator (LOL)" 75 if (ELFLOL_type==1) write(*,*) "10 Localized orbital locator (LOL) defined by Tsirelson" 76 if (ELFLOL_type==2) write(*,*) "10 Localized orbital locator (LOL) defined by Lu, Tian" 77 write(*,"(a,i5)") " 100 User-defined real space function, iuserfunc=",iuserfunc 78! write(*,*) "12 Total electrostatic potential" 79 read(*,*) ifunctopo 80 if (ifunctopo==4) then 81 write(*,"(a,i10)") " Input orbital index, between 1 and",nmo 82 read(*,*) iorbsel 83 end if 84 if (ifunctopo/=0) then 85 numcp=0 !Clean number 86 numpath=0 87 nple3n1path=0 88 numbassurf=0 89 CPtype=0 !Clean relationship 90 cp2surf=0 91 cp2ple3n1path=0 92 if (allocated(topopath)) deallocate(topopath) 93 if (allocated(bassurpath)) deallocate(bassurpath) 94 if (allocated(ple3n1path)) deallocate(ple3n1path) 95 write(*,*) "Note: All found CPs, paths, surfaces have been clean" 96 !Set special parameters for specific real space functions 97 if (ifunctopo==1.or.ifunctopo==4) then !High criteria for full analytical functions 98 gradconv=1D-7 99 dispconv=1D-8 100 else !User lower criteria for those functions with numerical gradient and Hessian 101 gradconv=1D-5 102 dispconv=1D-6 103 end if 104 if (ifunctopo==1) then 105 toposphrad=3D0 106 maxpathpt=451 !Default parameters for rho 107 pathstepsize=0.03D0 108 numsearchpt=2000 109 nsurfpathpercp=60 110 nsurfpt=100 111 surfpathstpsiz=0.03D0 112 else 113 toposphrad=3D0 114 if (ifunctopo==4) toposphrad=1.5D0 !CPs for orbital wavefunction is very close to nuclei, so use smaller radii 115 maxpathpt=901 116 pathstepsize=0.015D0 !Curvature is much larger than paths for rho, so smaller differentiate step must be used 117 numsearchpt=1000 118 if (ifunctopo==4) numsearchpt=10000 !Use large value since evaluation of orbital wavefunction is very fast, but hard to locate CPs 119 nsurfpathpercp=200 120 nsurfpt=100 121 surfpathstpsiz=0.008D0 122 end if 123 end if 124 else if (isel==-10) then 125 exit 126!-9 -9 -9 -9 -9 -9 -9 127 else if (isel==-9) then 128 write(*,*) "q = quit" 129 write(*,*) "Selection method: a? = Atom?, c? = Critical point ?" 130 write(*,*) "e.g. ""a1 c3"" returns the distance between atom1 and CP3" 131 write(*,*) " ""a4 a2"" returns the distance between atom4 and atom2" 132 write(*,*) " ""c6 a2 a5"" returns the angle of CP6-atom2-atom5" 133 write(*,*) " ""c2 c4 a3 c7"" returns the dihedral angle of CP2-CP4-atom3-CP7" 134 do while(.true.) 135 read(*,"(a)") c200 136 c200=adjustl(c200) 137 imeasure=0 138 do ichar=1,len_trim(c200) !imeasure=1/2/3: measure distance,angle,dihedral 139 if (c200(ichar:ichar)==','.or.c200(ichar:ichar)==' ') imeasure=imeasure+1 140 end do 141 if (c200(1:1)=='q') then 142 exit 143 else if (imeasure==1.or.imeasure==2.or.imeasure==3) then 144 if (imeasure==1) read(c200,*) ctmp1,ctmp2 !Read two terms 145 if (imeasure==2) read(c200,*) ctmp1,ctmp2,ctmp3 !Read three terms 146 if (imeasure==3) read(c200,*) ctmp1,ctmp2,ctmp3,ctmp4 !Read four terms 147 148 if (ctmp1(1:1)=='a') then 149 read(ctmp1(2:),*) iatm 150 tmpx1=a(iatm)%x 151 tmpy1=a(iatm)%y 152 tmpz1=a(iatm)%z 153 else if (ctmp1(1:1)=='c') then 154 read(ctmp1(2:),*) icp 155 tmpx1=CPpos(1,icp) 156 tmpy1=CPpos(2,icp) 157 tmpz1=CPpos(3,icp) 158 end if 159 if (ctmp2(1:1)=='a') then 160 read(ctmp2(2:),*) iatm 161 tmpx2=a(iatm)%x 162 tmpy2=a(iatm)%y 163 tmpz2=a(iatm)%z 164 else if (ctmp2(1:1)=='c') then 165 read(ctmp2(2:),*) icp 166 tmpx2=CPpos(1,icp) 167 tmpy2=CPpos(2,icp) 168 tmpz2=CPpos(3,icp) 169 end if 170 if (imeasure==1) write(*,"(' The distance is',f12.6,' Bohr (',f12.6 ' Angstrom)')") & 171 dsqrt((tmpx1-tmpx2)**2+(tmpy1-tmpy2)**2+(tmpz1-tmpz2)**2),dsqrt((tmpx1-tmpx2)**2+(tmpy1-tmpy2)**2+(tmpz1-tmpz2)**2)*b2a 172 173 if (imeasure==2.or.imeasure==3) then !Analyze one more term, then print angle 174 if (ctmp3(1:1)=='a') then 175 read(ctmp3(2:),*) iatm 176 tmpx3=a(iatm)%x 177 tmpy3=a(iatm)%y 178 tmpz3=a(iatm)%z 179 else if (ctmp3(1:1)=='c') then 180 read(ctmp3(2:),*) icp 181 tmpx3=CPpos(1,icp) 182 tmpy3=CPpos(2,icp) 183 tmpz3=CPpos(3,icp) 184 end if 185 end if 186 if (imeasure==2) write(*,"(' The angle is',f12.6,' degree')") xyz2angle(tmpx1,tmpy1,tmpz1,tmpx2,tmpy2,tmpz2,tmpx3,tmpy3,tmpz3) 187 188 if (imeasure==3) then !Analyze one more term, then print dihedral angle 189 if (ctmp4(1:1)=='a') then 190 read(ctmp4(2:),*) iatm 191 tmpx4=a(iatm)%x 192 tmpy4=a(iatm)%y 193 tmpz4=a(iatm)%z 194 else if (ctmp4(1:1)=='c') then 195 read(ctmp4(2:),*) icp 196 tmpx4=CPpos(1,icp) 197 tmpy4=CPpos(2,icp) 198 tmpz4=CPpos(3,icp) 199 end if 200 write(*,"(' The dihedral angle is',f12.6,' degree')") xyz2dih(tmpx1,tmpy1,tmpz1,tmpx2,tmpy2,tmpz2,tmpx3,tmpy3,tmpz3,tmpx4,tmpy4,tmpz4) 201 end if 202 else 203 write(*,*) "Input error" 204 end if 205 end do 206 207!-5 -5 -5 -5 -5 -5 -5 208!-5 -5 -5 -5 -5 -5 -5 209!-5 -5 -5 -5 -5 -5 -5 210 else if (isel==-5) then 211 do while(.true.) 212 write(*,*) 213 write(*,*) " ======== Process paths ========" 214 write(*,*) "0 Return" 215 write(*,*) "1 Print summary of paths" 216 write(*,*) "2 Print detail of a path" 217 write(*,*) "3 Delete some paths" 218 write(*,*) "4 Save points of all paths to paths.txt in current folder" 219 write(*,*) "5 Load paths from an external file" 220 write(*,*) "6 Export paths as paths.pdb file in current folder" 221 write(*,*) "7 Calculate specific real space funtion along the path" 222 read(*,*) isel2 223 224 if (isel2==0) then 225 exit 226 227 else if (isel2==1) then 228 if (numpath>0) then 229 do i=1,numpath 230 call path_cp(i,icp1,icp2,ipathtype) 231 if (icp1==0) then 232 icp1text=" Unknown " 233 else 234 write(icp1text,"(i5,1x,a)") icp1,CPtyp2lab(CPtype(icp1)) 235 end if 236 if (icp2==0) then 237 icp2text=" Unknown " 238 else 239 write(icp2text,"(i5,1x,a)") icp2,CPtyp2lab(CPtype(icp2)) 240 end if 241 write(*,"('#',i5,5x,'CP:',a,' --->',' CP:',a,' Length:',f9.5)") i,icp1text,icp2text,(pathnumpt(i)-1)*pathstepsize 242 end do 243 else 244 write(*,*) "No paths have been found" 245 end if 246 247 else if (isel2==2) then 248 write(*,*) "Input the index of the path" 249 read(*,*) ipath 250 if (ipath>numpath.or.ipath<=0) then 251 write(*,*) "Invalid index" 252 else 253 write(*,"(a,i6,a,f10.5,a,i5)") "Path:",ipath," Length:",(pathnumpt(ipath)-1)*pathstepsize," Total points:",pathnumpt(ipath) 254 write(*,"('From',3f18.12,/,'to ',3f18.12)") topopath(:,1,ipath),topopath(:,pathnumpt(ipath),ipath) 255 write(*,*) "The X/Y/Z coordinate (Bohr) and length of points in the path:" 256 do ipt=1,pathnumpt(ipath) 257 write(*,"(i6,3f16.10,f9.4)") ipt,topopath(:,ipt,ipath),(ipt-1)*pathstepsize 258 end do 259 end if 260 261 else if (isel2==3) then 262 write(*,*) "Input the index range that will be deleted, e.g. 3,10" 263 write(*,*) "Note: Input 0,0 can delete all paths" 264 read(*,*) idelstart,idelend 265 266 if (idelstart==0.and.idelend==0) then 267 numpath=0 268 if (allocated(topopath)) deallocate(topopath) 269 else if ( (idelstart>idelend).or.idelend>numpath.or.idelstart<=0) then 270 write(*,*) "Invalid input" 271 else 272 numdel=idelend-idelstart+1 273 nafter=numpath-idelend 274 topopath(:,:,idelstart:idelstart+nafter-1)=topopath(:,:,idelend+1:numpath) 275 pathnumpt(idelstart:idelstart+nafter-1)=pathnumpt(idelend+1:numpath) 276 numpath=numpath-numdel 277 write(*,"(' Now there are',i6,' paths left')") numpath 278 end if 279 280 else if (isel2==4) then 281 open(10,file="paths.txt",status="replace") 282 write(10,"(2i10)") numpath,maxpathpt 283 do ipath=1,numpath 284 write(10,"(/,'Path index:',i10)") ipath 285 write(10,"(i10)") pathnumpt(ipath) 286 do ipt=1,pathnumpt(ipath) 287 write(10,"(3E20.12)") topopath(:,ipt,ipath) 288 end do 289 end do 290 close(10) 291 write(*,*) "Done, path information have been saved to paths.txt in current folder" 292 write(*,"(a)") " Units are in Bohr. The first two numbers respectively denote the number of paths, and the maximal number of points in the paths" 293 294 else if (isel2==5) then 295 write(*,"(a)") " Note: The format of the input file must be identical to the one outputted by option 4" 296 if (numpath>0) write(*,*) "Note: After loading the file, all current paths will be clean" 297 write(*,*) "Input filename, e.g. c:\paths.txt" 298 read(*,*) c200 299 inquire(file=c200,exist=alive) 300 if (alive.eqv..false.) then 301 write(*,*) "File not found" 302 else 303 open(10,file=c200,status="old") 304 read(10,*) numpath,maxpathpt 305 if (allocated(topopath)) deallocate(topopath) 306 allocate(topopath(3,maxpathpt,numpath)) 307 do ipath=1,numpath 308 read(10,*) 309 read(10,*) 310 read(10,*) pathnumpt(ipath) 311 do ipt=1,pathnumpt(ipath) 312 read(10,*) topopath(:,ipt,ipath) 313 end do 314 end do 315 close(10) 316 write(*,*) "Done, path information have been recovered from the file" 317 end if 318 319 else if (isel2==6) then 320 open(10,file="paths.pdb",status="replace") 321 itmp=0 322 do ipath=1,numpath 323 do ipt=1,pathnumpt(ipath) 324 itmp=itmp+1 325 write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,2f6.2,10x,a2)") "HETATM",itmp,' '//"C "//' ',"PTH",'A',ipath,topopath(:,ipt,ipath)*b2a,1.0,0.0,"C " 326 end do 327 write(10,"('TER')") 328 end do 329 close(10) 330 write(*,*) "Done, path information have been saved to paths.pdb in current folder" 331 332 else if (isel2==7) then 333 do while(.true.) 334 write(*,*) 335 write(*,*) "Input index of a path, e.g. 3" 336 write(*,*) "Input ""q"" can return" 337 write(*,"(a)") " Hint: If input index of two paths (e.g. 6,7) emitted from the same (3,-1) CP, & 338 then the real space function along the combined paths will be outputted" 339 if (allocated(MOsym)) write(*,"(a)") " Hint: You can input ""s"" to choose which irreducible representations will be taken into account in the real space function evaluation" 340 read(*,"(a)") c200 341 342 if (index(c200,'q')/=0) then 343 exit 344 else if (index(c200,'s')/=0) then 345 call SelMO_IRREP 346 cycle 347 end if 348 349 itwopath=0 350 if (index(c200,',')==0) then !Only one path 351 read(c200,*) ipath 352 else !Two paths 353 itwopath=1 354 read(c200,*) ipath,jpath 355 end if 356 if (ipath>numpath.or.ipath<=0 .or. (itwopath==1.and.(jpath>numpath.or.jpath<=0)) ) then 357 write(*,*) "Error: Invalid index" 358 cycle 359 end if 360 if (itwopath==1.and.( topopath(1,1,ipath)/=topopath(1,1,jpath) .or. topopath(2,1,ipath)/=topopath(2,1,jpath) .or. topopath(3,1,ipath)/=topopath(3,1,jpath) )) then 361 write(*,*) "Error: The two paths are not emitted from the same (3,-1) critical point!" 362 cycle 363 end if 364 write(*,*) "Select the real space function to be calculated along the path" 365 call selfunc_interface(iselfunc) 366 !Store the values along the path into curvex and curvey, show them as text and curve map 367 npointcurve=pathnumpt(ipath) 368 if (itwopath==1) npointcurve=pathnumpt(ipath)+pathnumpt(jpath)-1 369 if (allocated(curvex)) deallocate(curvex) 370 if (allocated(curvey)) deallocate(curvey) 371 allocate(curvex(npointcurve),curvey(npointcurve)) 372 if (itwopath==0) then 373 do ipt=1,pathnumpt(ipath) 374 curvex(ipt)=(ipt-1)*pathstepsize 375 curvey(ipt)=calcfuncall(iselfunc,topopath(1,ipt,ipath),topopath(2,ipt,ipath),topopath(3,ipt,ipath)) 376 end do 377 else if (itwopath==1) then 378 ipttmp=0 379 do ipt=pathnumpt(ipath),1,-1 380 ipttmp=ipttmp+1 381 curvex(ipttmp)=(ipttmp-1)*pathstepsize 382 curvey(ipttmp)=calcfuncall(iselfunc,topopath(1,ipt,ipath),topopath(2,ipt,ipath),topopath(3,ipt,ipath)) 383 end do 384 do jpt=2,pathnumpt(jpath) !The first point of jpath is (3,-1), which has been included in ipath above 385 ipttmp=ipttmp+1 386 curvex(ipttmp)=(ipttmp-1)*pathstepsize 387 curvey(ipttmp)=calcfuncall(iselfunc,topopath(1,jpt,jpath),topopath(2,jpt,jpath),topopath(3,jpt,jpath)) 388 end do 389 end if 390 write(*,*) " Index X/Y/Z Coordinate (Bohr) Dist. Value" 391 if (itwopath==0) then 392 do ipt=1,pathnumpt(ipath) 393 write(*,"(i6,3f14.8,f9.4,E16.8)") ipt,topopath(:,ipt,ipath),curvex(ipt),curvey(ipt) 394 end do 395 icurve_vertlinex=0 396 else if (itwopath==1) then 397 ipttmp=0 398 do ipt=pathnumpt(ipath),1,-1 399 ipttmp=ipttmp+1 400 write(*,"(i6,3f14.8,f9.4,E16.8)") ipttmp,topopath(:,ipt,ipath),curvex(ipttmp),curvey(ipttmp) 401 end do 402 ibcptmp=ipttmp 403 do jpt=2,pathnumpt(jpath) !The first point of jpath is (3,-1), which has been included in ipath above 404 ipttmp=ipttmp+1 405 write(*,"(i6,3f14.8,f9.4,E16.8)") ipttmp,topopath(:,jpt,jpath),curvex(ipttmp),curvey(ipttmp) 406 end do 407 write(*,"(' Note: The point',i5,' corresponds to (3,-1) critical point')") ibcptmp 408 icurve_vertlinex=1 !Draw a vertical line to highlight BCP point 409 curve_vertlinex=curvex(ibcptmp) 410 write(*,*) "The dash line corresponds to the position of BCP" 411 end if 412 !Show the result as curve map 413 if (minval(curvey)>=0) then 414 curveymin=0 415 curveymax=1.1D0*maxval(curvey) 416 steplaby=curveymax/11 417 else 418 exty=(maxval(curvey)-minval(curvey))/10 419 curveymin=minval(curvey)-exty 420 curveymax=maxval(curvey)+exty 421 steplaby=exty 422 end if 423 steplabx=0.5D0 424 ilog10y=0 425 do while(.true.) 426 write(*,*) 427 write(*,*) "0 Return" 428 write(*,*) "1 Plot the graph again" 429 write(*,*) "2 Save the graph in current folder" 430 if (ilog10y==0) then 431 write(*,"(' 3 Set range of Y-axis, current:',f13.5,' to',f13.5,' Step:',f12.5)") curveymin,curveymax,steplaby 432 write(*,*) "4 Switch Y-axis type, current: linear scaling" 433 else if (ilog10y==1) then 434 write(*,"(' 3 Set range of Y-axis, current:',1PE14.6,' to',1PE14.6)") 10**curveymin,10**curveymax 435 write(*,*) "4 Switch Y-axis type, current: logarithmic scaling" 436 end if 437 write(*,*) "5 Export the data of the path shown above as pathvalue.txt in current folder" 438 read(*,*) iselcurve 439 440 if (iselcurve==0) then 441 exit 442 else if (iselcurve==1) then 443 else if (iselcurve==2) then 444 else if (iselcurve==3) then 445 if (ilog10y==0) then 446 write(*,*) "Input lower and upper limits and step between labels, e.g. 0,1.5,0.2" 447 read(*,*) curveymin,curveymax,steplaby 448 else if (ilog10y==1) then 449 write(*,*) "Input minimum and maximum value of Y axis e.g. -1,5 means from 10^-1 to 10^5" 450 read(*,*) curveymin,curveymax 451 end if 452 else if (iselcurve==4) then 453 if (ilog10y==0) then 454 ilog10y=1 455 curveymin=log10(curveymin) 456 curveymax=log10(curveymax) 457 else 458 ilog10y=0 459 curveymin=10**curveymin 460 curveymax=10**curveymax 461 end if 462 else if (iselcurve==5) then 463 open(10,file="pathvalue.txt",status="replace") 464 if (itwopath==0) then 465 do ipt=1,pathnumpt(ipath) 466 write(10,"(i6,3f14.8,f9.4,E16.8)") ipt,topopath(:,ipt,ipath),curvex(ipt),curvey(ipt) 467 end do 468 else if (itwopath==1) then 469 ipttmp=0 470 do ipt=pathnumpt(ipath),1,-1 471 ipttmp=ipttmp+1 472 write(10,"(i6,3f14.8,f9.4,E16.8)") ipttmp,topopath(:,ipt,ipath),curvex(ipttmp),curvey(ipttmp) 473 end do 474 do jpt=2,pathnumpt(jpath) !The first point of jpath is (3,-1), which has been included in ipath above 475 ipttmp=ipttmp+1 476 write(10,"(i6,3f14.8,f9.4,E16.8)") ipttmp,topopath(:,jpt,jpath),curvex(ipttmp),curvey(ipttmp) 477 end do 478 end if 479 close(10) 480 write(*,*) "Done! The data has been outputted to pathvalue.txt in current folder" 481 end if 482 end do 483 icurve_vertlinex=0 484 deallocate(curvex,curvey) 485 end do 486 end if 487 end do 488 489!-4 -4 -4 -4 -4 -4 -4 490!-4 -4 -4 -4 -4 -4 -4 491!-4 -4 -4 -4 -4 -4 -4 492 else if (isel==-4) then 493 do while(.true.) 494 write(*,*) " ============ Modify or export found CPs ============" 495 write(*,*) "-1 Print summary of CPs (in Angstrom)" 496 write(*,*) "0 Return" 497 write(*,*) "1 Print summary of CPs (in Bohr)" 498 write(*,*) "2 Delete some CPs" 499 write(*,*) "3 Add a CP artificially" 500 write(*,*) "4 Save CPs to CPs.txt in current folder" 501 write(*,*) "5 Load CPs from a file" 502 write(*,*) "6 Export CPs as CPs.pdb file in current folder" 503 read(*,*) isel2 504 505 if (isel2==0) then 506 exit 507 else if (isel2==1.or.isel2==-1) then 508 if (numcp>0) then 509 write(*,*) "Summary of found CPs:" 510 write(*,*) " Index Coordinate Type" 511 do icp=1,numcp 512 if (isel2==1) write(*,"(i6,3f16.9,3x,a)") icp,CPpos(:,icp),CPtyp2lab(CPtype(icp)) 513 if (isel2==-1) write(*,"(i6,3f16.9,3x,a)") icp,CPpos(:,icp)*b2a,CPtyp2lab(CPtype(icp)) 514 end do 515 else 516 write(*,*) "No CPs have been found" 517 end if 518 else if (isel2==2) then 519 if (numbassurf>0) then 520 write(*,"(a)") "Warning: If one or more CPs are deleted, all interbasin surfaces will be lost, continue? 0/1=No/Yes" 521 read(*,*) ifok 522 if (ifok==0) then 523 cycle 524 else 525 numbassurf=0 526 cp2surf=0 527 deallocate(bassurpath) 528 end if 529 end if 530 write(*,*) "Input the index range that will be deleted, e.g. 3,10" 531 write(*,*) "Note 1: Input 0,0 can delete all CPs" 532 write(*,*) "Note 2: If you want to delete CPs with too low electron density, select -1" 533 read(*,*) idelstart,idelend 534 if (idelstart==0.and.idelend==0) then 535 numcp=0 536 else if ( (idelstart>idelend).or.idelend>numcp.or.idelstart<=0) then 537 write(*,*) "Invalid input" 538 else 539 numdel=idelend-idelstart+1 540 nafter=numcp-idelend 541 CPpos(:,idelstart:idelstart+nafter-1)=CPpos(:,idelend+1:numcp) 542 CPtype(idelstart:idelstart+nafter-1)=CPtype(idelend+1:numcp) 543 numcp=numcp-numdel 544 end if 545 write(*,"(' Now there are',i6,' CPs left')") numcp 546 else if (isel2==3) then 547 numcp=numcp+1 548 write(*,*) "Input coordinate, e.g. 1.2,0.0,3.44" 549 read(*,*) CPpos(1,numcp),CPpos(2,numcp),CPpos(3,numcp) 550 write(*,*) "Input CP type, 1=(3,-3) 2=(3,-1) 3=(3,+1) 4=(3,+3)" 551 read(*,*) CPtype(numcp) 552 else if (isel2==4) then 553 open(10,file="CPs.txt",status="replace") 554 write(10,"(i6)") numcp 555 do icp=1,numcp 556 write(10,"(i6,3f12.6,3x,i4)") icp,CPpos(:,icp),CPtype(icp) 557 end do 558 close(10) 559 write(*,*) "Done, CP information have been saved to CPs.txt in current folder" 560 write(*,*) "Note: The last column is CP type, 1=(3,-3) 2=(3,-1) 3=(3,+1) 4=(3,+3)" 561 else if (isel2==5) then 562 write(*,*) "Input filename, e.g. C:\ltwd\CPs.txt" 563 write(*,"(a)") " (The format of the file must be identical to the one outputted by option 4)" 564 if (numcp>0) write(*,*) "Note: After loading the file, all found CPs will be clean" 565 read(*,*) c200 566 inquire(file=c200,exist=alive) 567 if (alive.eqv..false.) then 568 write(*,*) "File not found" 569 else 570 open(10,file=c200,status="old") 571 read(10,*) numcp 572 do icp=1,numcp 573 read(10,*) nouse,CPpos(:,icp),CPtype(icp) 574 end do 575 close(10) 576 write(*,*) "Done, CP information have been recovered from the file" 577 end if 578 else if (isel2==6) then 579 open(10,file="CPs.pdb",status="replace") 580 do icp=1,numcp 581 if (CPtype(icp)==1) write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,2f6.2,10x,a2)") "HETATM",icp,' '//"C "//' ',"CPS",'A',1,CPpos(:,icp)*b2a,1.0,0.0,"C " 582 if (CPtype(icp)==2) write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,2f6.2,10x,a2)") "HETATM",icp,' '//"N "//' ',"CPS",'A',1,CPpos(:,icp)*b2a,1.0,0.0,"N " 583 if (CPtype(icp)==3) write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,2f6.2,10x,a2)") "HETATM",icp,' '//"O "//' ',"CPS",'A',1,CPpos(:,icp)*b2a,1.0,0.0,"O " 584 if (CPtype(icp)==4) write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,2f6.2,10x,a2)") "HETATM",icp,' '//"F "//' ',"CPS",'A',1,CPpos(:,icp)*b2a,1.0,0.0,"F " 585 end do 586 write(*,*) "Done, CP information have been saved to CPs.pdb in current folder" 587 write(*,*) "Note: Element C/N/O/F correspond to (3,-3)/(3,-1)/(3,+1)/(3,+3) respectively" 588 close(10) 589 end if 590 write(*,*) 591 end do 592!-3 -3 -3 -3 -3 -3 -3 593 else if (isel==-3) then 594 do while(.true.) 595 write(*,*) "============ Set interbasin surface generating parameters ============" 596 write(*,"(a)") " 0 Return" 597 write(*,"(a,i6)") " 1 Number of paths in each interbasin surface",nsurfpathpercp 598 write(*,"(a,i6)") " 2 Number of points in each interbasin surface path, current:",nsurfpt 599 write(*,"(a,f8.4)") " 3 Stepsize, current:",surfpathstpsiz 600 read(*,*) isel2 601 602 if (isel2==0) then 603 exit 604 else 605 if (allocated(bassurpath)) then 606 write(*,*) "If the parameter is changed, all already generated surfaces will be clean, OK?" 607 write(*,*) "0/1=No/Yes" 608 read(*,*) ifok 609 if (ifok==1) then 610 deallocate(bassurpath) 611 numbassurf=0 612 CP2surf=0 613 end if 614 end if 615 if (.not.allocated(bassurpath)) then 616 write(*,*) "Input the value" 617 if (isel2==1) read(*,*) nsurfpathpercp 618 if (isel2==2) read(*,*) nsurfpt 619 if (isel2==3) read(*,*) surfpathstpsiz 620 end if 621 end if 622 end do 623!-2 -2 -2 -2 -2 -2 -2 624 else if (isel==-2) then 625 do while(.true.) 626 write(*,*) " ============ Set path searching parameters ============" 627 write(*,"(a)") " 0 Return" 628 write(*,"(a,i4)") " 1 Maximal number of points of a path, current:",maxpathpt 629 write(*,"(a,f8.4)") " 2 Stepsize, current:",pathstepsize 630 write(*,"(a,f8.4)") " 3 Stop generation if distance to any CP is smaller than:",discritpathfin 631! write(*,"(a,i4)") " 4 Maximal number of bisection, current:",npathtry 632 read(*,*) isel2 633 634 if (isel2==0) then 635 exit 636 else if (isel2==1.or.isel2==2) then 637 iok=1 638 if (numpath>0) then 639 write(*,*) "Warning: All found paths will be clean, OK? 1=Yes 2=No" 640 read(*,*) ifok 641 end if 642 if (iok==1) then 643 write(*,*) "Input a value" 644 if (isel2==1) read(*,*) maxpathpt 645 if (isel2==2) read(*,*) pathstepsize 646 numpath=0 647 if (allocated(topopath)) deallocate(topopath) 648 end if 649 else if (isel2==3) then 650 write(*,*) "Input a value" 651 read(*,*) discritpathfin 652 else if (isel2==4) then 653 write(*,*) "Input an integer, must >=2" 654 read(*,*) npathtry 655 end if 656 end do 657!-1 -1 -1 -1 -1 -1 -1 658 else if (isel==-1) then 659 do while(.true.) 660 write(*,*)" ============ Set CP searching parameters ============" 661 write(*,"(a)") " -1 Set to default" 662 write(*,"(a)") " 0 Return" 663 write(*,"(a,i5)") " 1 Set maximal iterations:",topomaxcyc 664 write(*,"(a,f12.6)") " 2 Set scale factor for stepsize:",CPstepscale 665 write(*,"(a,1PE12.5)") " 3 Criteria for gradient-norm convergence:",gradconv 666 write(*,"(a,1PE12.5)") " 4 Criteria for displacement convergence:",dispconv 667 write(*,"(a,f12.6)") " 5 Minimal distance between CPs:",minicpdis 668 write(*,"(a,f8.2)") " 6 Skip search if distance between atoms is longer than the sum of their vdW radius multiplied by:",vdwsumcrit 669 if (ishowsearchlevel==0) write(*,"(a)") " 7 If print details of CP searching procedure: No" 670 if (ishowsearchlevel==1) write(*,"(a)") " 7 If print details of CP searching procedure: Some detail" 671 if (ishowsearchlevel==2) write(*,"(a)") " 7 If print details of CP searching procedure: All detail" 672 write(*,"(a,1PE15.8)") " 8 Criteria for determining if Hessian matrix is singular:",singularcrit 673 if (CPsearchlow==CPsearchhigh) then 674 write(*,*) "9 Select rule for reserving CPs, current: reserve All CPs" 675 else 676 write(*,"(a,1PE12.4,a,1PE12.4)") " 9 Select rule for reserving CPs, current: between",CPsearchlow,' and',CPsearchhigh 677 end if 678 read(*,*) isel2 679 680 if (isel2==-1) then 681 topomaxcyc=120 682 CPstepscale=1D0 683 gradonv=1D-7 684 dispconv=1D-9 685 minicpdis=0.03D0 686 vdwsumcrit=1.2D0 687 ishowsearchlevel=0 688 else if (isel2==0) then 689 exit 690 else if (isel2==1) then 691 write(*,*) "Input an integer" 692 read(*,*) topomaxcyc 693 else if (isel2==2) then 694 write(*,*) "Input a value, default is 1.0" 695 read(*,*) CPstepscale 696 else if (isel2==3) then 697 write(*,*) "Input a value" 698 read(*,*) gradconv 699 else if (isel2==4) then 700 write(*,*) "Input a value" 701 read(*,*) dispconv 702 else if (isel2==5) then 703 write(*,*) "Input a value" 704 read(*,*) minicpdis 705 else if (isel2==6) then 706 write(*,*) "Input a value" 707 read(*,*) vdwsumcrit 708 else if (isel2==7) then 709 write(*,*) "0 Don't print details" 710 write(*,*) "1 Print some details" 711 write(*,*) "2 Print all details" 712 read(*,*) ishowsearchlevel 713 if ( nthreads >1) write(*,*) "Warning: The printed details may be messed up since parallel mode is enabled!" 714 else if (isel2==8) then 715 write(*,"(a)") "Input a value, if absolute value of determinant of Hessiant matrix is lower than this value, then it will be regarded as singular, e.g. 1D-21" 716 read(*,*) singularcrit 717 else if (isel2==9) then 718 write(*,"(a)") " Input lower and upper limits. For example, if you input 0.05,0.22, & 719 then during CP searching, only when the real space function of a new CP is between 0.05 and 0.22 then it will be reserved" 720 write(*,*) "Note: If the two values are identical, all CPs will be reserved (default case)" 721 read(*,*) CPsearchlow,CPsearchhigh 722 end if 723 end do 724!0000000000000000 725 else if (isel==0) then 726 if (numpath>0) then 727 write(*,*) "Summary of found paths:" 728 do i=1,numpath 729 call path_cp(i,icp1,icp2,ipathtype) 730 if (icp1==0) then 731 icp1text=" Unknown " 732 else 733 write(icp1text,"(i5,1x,a)") icp1,CPtyp2lab(CPtype(icp1)) 734 end if 735 if (icp2==0) then 736 icp2text=" Unknown " 737 else 738 write(icp2text,"(i5,1x,a)") icp2,CPtyp2lab(CPtype(icp2)) 739 end if 740 write(*,"('#',i5,5x,'CP:',a,' --->',' CP:',a,' Length:',f9.5)") i,icp1text,icp2text,(pathnumpt(i)-1)*pathstepsize 741 end do 742 else 743 write(*,*) "No paths have been found" 744 end if 745 write(*,*) 746 if (numcp>0) then 747 write(*,*) "Summary of found CPs:" 748 write(*,*) " Index XYZ Coordinate (Bohr) Type" 749 do icp=1,numcp 750 icptype=CPtype(icp) 751 if (ifunctopo==1.and.icptype==1) then 752 do iatm=1,ncenter 753 disttmp2=(CPpos(1,icp)-a(iatm)%x)**2+(CPpos(2,icp)-a(iatm)%y)**2+(CPpos(3,icp)-a(iatm)%z)**2 754 if (disttmp2<0.01D0) then 755 write(*,"(i6,3f16.9,3x,a,' Nuc:',i5,'(',a')')") icp,CPpos(:,icp),CPtyp2lab(icptype),iatm,a(iatm)%name 756 exit 757 end if 758 end do 759 if (iatm==ncenter+1) write(*,"(i6,3f16.9,3x,a,' Nuc: Unknown')") icp,CPpos(:,icp),CPtyp2lab(icptype) 760 else 761 write(*,"(i6,3f16.9,3x,a)") icp,CPpos(:,icp),CPtyp2lab(icptype) 762 end if 763 end do 764 else 765 write(*,*) "No CPs have been found" 766 end if 767 NumCPtype1=count(CPtype(1:numcp)==1) 768 NumCPtype2=count(CPtype(1:numcp)==2) 769 NumCPtype3=count(CPtype(1:numcp)==3) 770 NumCPtype4=count(CPtype(1:numcp)==4) 771 write(*,*) "The number of critical points of each type:" 772 write(*,"(' (3,-3):',i6,', (3,-1):',i6,', (3,+1):',i6,', (3,+3):',i6)") NumCPtype1,NumCPtype2,NumCPtype3,NumCPtype4 773 774 itestPH=NumCPtype1-NumCPtype2+NumCPtype3-NumCPtype4 !Poincare-Hopf relationship 775 write(*,"(' Poincare-Hopf relationship verification:',i5,' -',i5,' +',i5,' -',i5,' =',i4)") NumCPtype1,NumCPtype2,NumCPtype3,NumCPtype4,itestPH 776 if (itestPH/=1) write(*,*) "Warning: Poincare-Hopf relationship is not satisfied, some CPs may be missing" 777 if (itestPH==1) write(*,*) "Fine, Poincare-Hopf relationship is satisfied, all CPs may have been found" 778 779 if (numbassurf>0) write(*,"(' The number of generated interbasin surfaces:',i8)") numbassurf 780 if (numpath>0) idrawmol=0 !Avoid atom and bond covered paths 781!111111111111111111111 782 else if (isel==1) then 783 numcpold=numcp 784 write(*,*) "Input X,Y,Z of starting point (in bohr, e.g. 2.0,3.1,-0.5)" 785 write(*,"(a)") " You can also input two atomic indices (e.g. 4,5), then midpoint of corresponding two atoms will be taken as starting point" 786 read(*,"(a)") c200 787 read(c200,*,iostat=ierror) x,y,z 788 if (ierror/=0) then 789 read(c200,*) iatm,jatm 790 x=(a(iatm)%x+a(jatm)%x)/2 791 y=(a(iatm)%y+a(jatm)%y)/2 792 z=(a(iatm)%z+a(jatm)%z)/2 793 end if 794 call findcp(x,y,z,ifunctopo,0) 795 if (numcp==numcpold) then 796 write(*,*) "No new critical point was found" 797 else 798 write(*,"(' Find',i5,' new critical points')") numcp-numcpold 799 end if 800!222222222222222222222 801 else if (isel==2) then 802 numcpold=numcp 803 do iatm=1,ncenter 804 write(*,"('#',i5,' /',i5,a,i5,'(',a,')')") iatm,ncenter,": Trying from nuclear position of ",iatm,a(iatm)%name 805 if (a(iatm)%index>10) then 806 call findcp(a(iatm)%x,a(iatm)%y,a(iatm)%z,ifunctopo,1) !For heavy atoms, use lower criteria, because the cusp of electron density is sharp so hard to locate 807 else 808 call findcp(a(iatm)%x,a(iatm)%y,a(iatm)%z,ifunctopo,0) 809 end if 810 end do 811 if ((numcp-numcpold)/=0) then 812 call sortCP(numcpold+1) 813 write(*,*) " ==== Summary ====" 814 write(*,*) " Index Coordinate Type" 815 do icp=numcpold+1,numcp 816 write(*,"(i6,3f12.6,3x,a)") icp,CPpos(:,icp),CPtyp2lab(CPtype(icp)) 817 end do 818 end if 819 write(*,"(' Totally find',i6,' new critical points')") numcp-numcpold 820 if (ifunctopo==1.and.count(CPtype(1:numcp)==1)<ncenter) write(*,*) "Warning: Some (3,-3) may missing, try to search again with different parameters" 821!333333333333333333333 822 else if (isel==3) then 823 numcpold=numcp 824 itime=0 825 ntime=0 826 do iatm=1,ncenter !Test how many iterations will be done 827 do jatm=iatm+1,ncenter 828 if ( distmat(iatm,jatm) <= vdwsumcrit*(vdwr(a(iatm)%index)+vdwr(a(jatm)%index)) ) ntime=ntime+1 829 end do 830 end do 831nthreads=getNThreads() 832!$OMP PARALLEL DO SHARED(itime) PRIVATE(iatm,jatm) schedule(dynamic) NUM_THREADS(nthreads) 833 do iatm=1,ncenter 834 do jatm=iatm+1,ncenter 835 if ( distmat(iatm,jatm) > vdwsumcrit*(vdwr(a(iatm)%index)+vdwr(a(jatm)%index)) ) cycle 836!$OMP CRITICAL 837 itime=itime+1 838!$OMP end CRITICAL 839 write(*,"('#',i5,' /',i5,a,i5,'(',a,')',a,i5,'(',a,')')") & 840 itime,ntime,": Trying from midpoint between ",iatm,a(iatm)%name," and",jatm,a(jatm)%name 841 call findcp( (a(iatm)%x+a(jatm)%x)/2D0,(a(iatm)%y+a(jatm)%y)/2D0,(a(iatm)%z+a(jatm)%z)/2D0, ifunctopo,0) 842 end do 843 end do 844!$OMP END PARALLEL DO 845 if ((numcp-numcpold)/=0) then 846 call sortCP(numcpold+1) 847 write(*,*) " ==== Summary ====" 848 write(*,*) " Index Coordinate Type" 849 do icp=numcpold+1,numcp 850 write(*,"(i6,3f12.6,3x,a)") icp,CPpos(:,icp),CPtyp2lab(CPtype(icp)) 851 end do 852 end if 853 write(*,"(' Totally find',i6,' new critical points')") numcp-numcpold 854!4444444444444444444 855 else if (isel==4) then 856 numcpold=numcp 857 itime=0 858 ntime=0 859 do iatm=1,ncenter !Test how many iterations will be done 860 do jatm=iatm+1,ncenter 861 if ( distmat(iatm,jatm) > vdwsumcrit*(vdwr(a(iatm)%index)+vdwr(a(jatm)%index)) ) cycle 862 do katm=jatm+1,ncenter 863 if ( distmat(katm,iatm) > vdwsumcrit*(vdwr(a(katm)%index)+vdwr(a(iatm)%index)).or.& 864 distmat(katm,jatm) > vdwsumcrit*(vdwr(a(katm)%index)+vdwr(a(jatm)%index)) ) cycle 865 ntime=ntime+1 866 end do 867 end do 868 end do 869nthreads=getNThreads() 870!$OMP PARALLEL DO SHARED(itime) PRIVATE(iatm,jatm,katm) schedule(dynamic) NUM_THREADS(nthreads) 871 do iatm=1,ncenter 872 do jatm=iatm+1,ncenter 873 if ( distmat(iatm,jatm) > vdwsumcrit*(vdwr(a(iatm)%index)+vdwr(a(jatm)%index)) ) cycle 874 do katm=jatm+1,ncenter 875 if ( distmat(katm,iatm) > vdwsumcrit*(vdwr(a(katm)%index)+vdwr(a(iatm)%index)).or.& 876 distmat(katm,jatm) > vdwsumcrit*(vdwr(a(katm)%index)+vdwr(a(jatm)%index)) ) cycle 877!$OMP CRITICAL 878 itime=itime+1 879!$OMP end CRITICAL 880 write(*,"('#',i5,' /',i5,a ,i5,'(',a,')' ,i5,'(',a,')' ,i5,'(',a,')' )") & 881 itime,ntime,": Trying from triangle center of ",iatm,a(iatm)%name,jatm,a(jatm)%name,katm,a(katm)%name 882 call findcp( (a(iatm)%x+a(jatm)%x+a(katm)%x)/3D0,(a(iatm)%y+a(jatm)%y+a(katm)%y)/3D0,(a(iatm)%z+a(jatm)%z+a(katm)%z)/3D0, ifunctopo,0) 883 end do 884 end do 885 end do 886!$OMP END PARALLEL DO 887 if ((numcp-numcpold)/=0) then 888 call sortCP(numcpold+1) 889 write(*,*) " ==== Summary ====" 890 write(*,*) " Index Coordinate Type" 891 do icp=numcpold+1,numcp 892 write(*,"(i6,3f12.6,3x,a)") icp,CPpos(:,icp),CPtyp2lab(CPtype(icp)) 893 end do 894 end if 895 write(*,"(' Totally find',i6,' new critical points')") numcp-numcpold 896!5555555555555555555 897 else if (isel==5) then 898 numcpold=numcp 899 itime=0 900 ntime=0 901 do iatm=1,ncenter !Test how many iterations will be done; ij,jk,kl,li,lj,ik 902 do jatm=iatm+1,ncenter 903 if ( distmat(iatm,jatm) > vdwsumcrit*(vdwr(a(iatm)%index)+vdwr(a(jatm)%index)) ) cycle 904 do katm=jatm+1,ncenter 905 if ( distmat(katm,jatm) > vdwsumcrit*(vdwr(a(katm)%index)+vdwr(a(jatm)%index)) ) cycle 906 do latm=katm+1,ncenter 907 if ( distmat(latm,katm) > vdwsumcrit*(vdwr(a(latm)%index)+vdwr(a(katm)%index)).or.& 908 distmat(latm,iatm) > vdwsumcrit*(vdwr(a(latm)%index)+vdwr(a(iatm)%index)).or.& 909 distmat(latm,jatm) > vdwsumcrit*(vdwr(a(latm)%index)+vdwr(a(jatm)%index)).or.& 910 distmat(iatm,katm) > vdwsumcrit*(vdwr(a(iatm)%index)+vdwr(a(katm)%index))) cycle 911 ntime=ntime+1 912 end do 913 end do 914 end do 915 end do 916nthreads=getNThreads() 917!$OMP PARALLEL DO SHARED(itime) PRIVATE(iatm,jatm,katm,latm) schedule(dynamic) NUM_THREADS(nthreads) 918 do iatm=1,ncenter !Test how many iterations will be done; ij,jk,kl,li,lj,ik 919 do jatm=iatm+1,ncenter 920 if ( distmat(iatm,jatm) > vdwsumcrit*(vdwr(a(iatm)%index)+vdwr(a(jatm)%index)) ) cycle 921 do katm=jatm+1,ncenter 922 if ( distmat(katm,jatm) > vdwsumcrit*(vdwr(a(katm)%index)+vdwr(a(jatm)%index)) ) cycle 923 do latm=katm+1,ncenter 924 if ( distmat(latm,katm) > vdwsumcrit*(vdwr(a(latm)%index)+vdwr(a(katm)%index)).or.& 925 distmat(latm,iatm) > vdwsumcrit*(vdwr(a(latm)%index)+vdwr(a(iatm)%index)).or.& 926 distmat(latm,jatm) > vdwsumcrit*(vdwr(a(latm)%index)+vdwr(a(jatm)%index)).or.& 927 distmat(iatm,katm) > vdwsumcrit*(vdwr(a(iatm)%index)+vdwr(a(katm)%index))) cycle 928!$OMP CRITICAL 929 itime=itime+1 930!$OMP end CRITICAL 931 write(*,"('#',i5,' /',i5,a ,i5,'(',a,')' ,i5,'(',a,')' ,i5,'(',a,')' ,i5,'(',a,')')") & 932 itime,ntime,": Trying from center of ",& 933 iatm,a(iatm)%name,jatm,a(jatm)%name,katm,a(katm)%name,latm,a(latm)%name 934 call findcp( (a(iatm)%x+a(jatm)%x+a(katm)%x+a(latm)%x)/4D0,& 935 (a(iatm)%y+a(jatm)%y+a(katm)%y+a(latm)%y)/4D0,(a(iatm)%z+a(jatm)%z+a(katm)%z+a(latm)%z)/4D0, ifunctopo,0) 936 end do 937 end do 938 end do 939 end do 940!$OMP END PARALLEL DO 941 if ((numcp-numcpold)/=0) then 942 call sortCP(numcpold+1) 943 write(*,*) " ==== Summary ====" 944 write(*,*) " Index Coordinate Type" 945 do icp=numcpold+1,numcp 946 write(*,"(i6,3f12.6,3x,a)") icp,CPpos(:,icp),CPtyp2lab(CPtype(icp)) 947 end do 948 end if 949 write(*,"(' Totally found',i6,' new critical points')") numcp-numcpold 950!6666666666666666666 951 else if (isel==6) then 952 do while(.true.) 953 write(*,*) " ============= Set the starting points =============" 954 write(*,"(' Center:',3f10.5,' Radius:',f6.2,' Points:',i8)") sphcenx,sphceny,sphcenz,toposphrad,numsearchpt 955 write(*,"(a)") " -9 Return" 956 write(*,"(a)") " -2 Start the search using some nuclei as sphere center in turn" 957 write(*,"(a)") " -1 Start the search using each nucleus as sphere center in turn" 958 write(*,"(a)") " 0 Start the search using the defined sphere center" 959 write(*,"(a)") " 1 Input coordinate of the sphere center" 960 write(*,"(a)") " 2 Set the sphere center at a nuclear position" 961 write(*,"(a)") " 3 Set the sphere center at midpoint between two atoms" 962 write(*,"(a)") " 4 Set the sphere center at triangle center of three atoms" 963 write(*,"(a)") " 5 Set the sphere center at a CP" 964 write(*,"(a)") " 6 Set the sphere center at midpoint between two CPs" 965 write(*,"(a)") " 10 Set the sphere radius" 966 write(*,"(a)") " 11 Set the number of points in the sphere" 967 read(*,*) isel2 968 969 if (isel2==-9) then 970 exit 971 else if (isel2==0.or.isel2==-1.or.isel2==-2) then 972 if (isel2==-2) then 973 write(*,*) "Input the indices of the atoms, e.g. 3,4,5,12, at most 1000 characters" 974 read(*,"(a)") c1000 975 call str2arr(c1000,nsearchcen,searchcenlist) 976 end if 977 if (isel2==0) nsearchcen=1 978 if (isel2==-1) nsearchcen=ncenter 979 !4.189=4/3*pi, this assess how many points need to be generated in the cube 980 !so that numsearchpt points could in the sphere 981 numcpold=numcp 982 numsearchpt_tmp=nint(8D0/4.189D0*numsearchpt) 983 allocate(randptx(numsearchpt_tmp),randpty(numsearchpt_tmp),randptz(numsearchpt_tmp)) 984 985 itime=0 986 ioutcount=0 987 do icenidx=1,nsearchcen 988 icen=icenidx !isel==0.or.isel==-1 989 if (isel2==-2) icen=searchcenlist(icenidx) 990 if (isel2==-1.or.isel2==-2) then !Cycle each atom center 991 sphcenx=a(icen)%x 992 sphceny=a(icen)%y 993 sphcenz=a(icen)%z 994 end if 995 CALL RANDOM_NUMBER(randptx) 996 CALL RANDOM_NUMBER(randpty) 997 CALL RANDOM_NUMBER(randptz) 998 randptx=randptx*2*toposphrad+(sphcenx-toposphrad) !Move distribution center of random point to sphere center 999 randpty=randpty*2*toposphrad+(sphceny-toposphrad) 1000 randptz=randptz*2*toposphrad+(sphcenz-toposphrad) 1001 randptx(1)=sphcenx !The first try point is set to sphere center, this is faciliate to locate CP at nuclei 1002 randpty(1)=sphceny 1003 randptz(1)=sphcenz 1004nthreads=getNThreads() 1005!$OMP PARALLEL DO SHARED(itime) PRIVATE(i) schedule(dynamic) NUM_THREADS(nthreads) 1006 do i=1,numsearchpt_tmp 1007 dispt_cen=dsqrt( (randptx(i)-sphcenx)**2+(randpty(i)-sphceny)**2+(randptz(i)-sphcenz)**2 ) 1008 if (dispt_cen>toposphrad) cycle 1009 call findcp(randptx(i),randpty(i),randptz(i),ifunctopo,0) 1010!$OMP CRITICAL 1011 itime=itime+1 1012 if (itime>ioutcount+99.or.ifunctopo==12) then 1013 write(*,"('#',i10,' /',i10)") itime,numsearchpt*nsearchcen 1014 ioutcount=ioutcount+100 1015 end if 1016!$OMP end CRITICAL 1017 end do 1018!$OMP END PARALLEL DO 1019 end do 1020 deallocate(randptx,randpty,randptz) 1021 1022 if ((numcp-numcpold)/=0) then 1023! call sortCP(numcpold+1) !Senseless here, because the guessing points occur randomly 1024 write(*,*) " ==== Summary ====" 1025 write(*,*) " Index Coordinate Type" 1026 do icp=numcpold+1,numcp 1027 write(*,"(i6,3f12.6,3x,a)") icp,CPpos(:,icp),CPtyp2lab(CPtype(icp)) 1028 end do 1029 end if 1030 write(*,"(' Totally find',i6,' new critical points')") numcp-numcpold 1031 write(*,*) 1032 else if (isel2==1) then 1033 write(*,*) "Input x,y,z e.g. 1.2,0.2,-0.44" 1034 read(*,*) sphcenx,sphceny,sphcenz 1035 else if (isel2==2) then 1036 write(*,*) "Input atom index" 1037 read(*,*) iatm 1038 if (iatm>ncenter.or.iatm<=0) then 1039 write(*,*) "Invalid input" 1040 else 1041 sphcenx=a(iatm)%x 1042 sphceny=a(iatm)%y 1043 sphcenz=a(iatm)%z 1044 end if 1045 else if (isel2==3) then 1046 write(*,*) "Input index of the two atoms, e.g. 3,7" 1047 read(*,*) iatm,jatm 1048 if ( iatm>ncenter.or.iatm<=0.or.jatm>ncenter.or.jatm<=0 ) then 1049 write(*,*) "Invalid input" 1050 else 1051 sphcenx=(a(iatm)%x+a(jatm)%x)/2D0 1052 sphceny=(a(iatm)%y+a(jatm)%y)/2D0 1053 sphcenz=(a(iatm)%z+a(jatm)%z)/2D0 1054 end if 1055 else if (isel2==4) then 1056 write(*,*) "Input index of the three atoms, e.g. 2,3,7" 1057 read(*,*) iatm,jatm,katm 1058 if ( iatm>ncenter.or.iatm<=0.or.jatm>ncenter.or.jatm<=0.or.katm>ncenter.or.katm<=0 ) then 1059 write(*,*) "Invalid input" 1060 else 1061 sphcenx=(a(iatm)%x+a(jatm)%x+a(katm)%x)/3D0 1062 sphceny=(a(iatm)%y+a(jatm)%y+a(katm)%y)/3D0 1063 sphcenz=(a(iatm)%z+a(jatm)%z+a(katm)%z)/3D0 1064 end if 1065 else if (isel2==5) then 1066 write(*,*) "Input CP index" 1067 read(*,*) icp 1068 if (icp>numcp.or.icp<=0) then 1069 write(*,*) "Invalid input" 1070 else 1071 sphcenx=CPpos(1,icp) 1072 sphceny=CPpos(2,icp) 1073 sphcenz=CPpos(3,icp) 1074 end if 1075 else if (isel2==6) then 1076 write(*,*) "Input index of the two CPs, e.g. 3,7" 1077 read(*,*) icp,jcp 1078 if ( icp>numcp.or.icp<=0.or.jcp>numcp.or.jcp<=0 ) then 1079 write(*,*) "Invalid input" 1080 else 1081 sphcenx=(CPpos(1,icp)+CPpos(1,jcp))/2D0 1082 sphceny=(CPpos(2,icp)+CPpos(2,jcp))/2D0 1083 sphcenz=(CPpos(3,icp)+CPpos(3,jcp))/2D0 1084 end if 1085 else if (isel2==10) then 1086 write(*,*) "Input a radius" 1087 read(*,*) toposphrad 1088 else if (isel2==11) then 1089 write(*,*) "Input a number" 1090 read(*,*) numsearchpt 1091 end if 1092 end do 1093!7777777777777777777 1094 else if (isel==7) then 1095 write(*,*) "Input the index of the CP that you are interested in, e.g. 3" 1096 write(*,"(a)") " Note 1: If input 0, then properties of all CPs will be outputted to CPprop.txt in current folder & 1097 (and if you feel the output speed is slow, you can input -1 to avoid outputting ESP, which is the most expensive one)" 1098 write(*,"(a)") " Note 2: If input CP index with ""d"" suffix, e.g. 17d, then property of this CP can be decomposed into orbital contribution" 1099 read(*,*) c200 1100 if (index(c200,"d")/=0) then 1101 read(c200(1:index(c200,"d")-1),*) indcp 1102 call decompptprop(CPpos(1,indcp),CPpos(2,indcp),CPpos(3,indcp)) 1103 else 1104 read(c200,*) indcp 1105 if (indcp==0.or.indcp==-1) then 1106 write(*,*) "Please wait..." 1107 iback=ishowptESP 1108 if (indcp==-1) ishowptESP=0 !Avoid outputting ESP 1109 open(10,file="CPprop.txt",status="replace") 1110 do icp=1,numcp 1111 write(*,"(' Outputting CP',i6,' /',i6)") icp,numcp 1112 write(10,"(' ================ CP',i6,', Type ',a,' ================')") icp,CPtyp2lab(CPtype(icp)) 1113 write(10,"(' Position (Bohr):',3f20.14)") CPpos(:,icp) 1114 call showptprop(CPpos(1,icp),CPpos(2,icp),CPpos(3,icp),ifunctopo,10) 1115 write(10,*) 1116 end do 1117 close(10) 1118 if (indcp==-1) ishowptESP=iback 1119 write(*,*) "Done! The results have been outputted to CPprop.txt in current folder" 1120 write(*,*) "Note: Unless otherwise specified, all units are in a.u." 1121 else if (indcp>0.and.indcp<=numcp) then 1122 write(*,*) "Note: Unless otherwise specified, all units are in a.u." 1123 write(*,"(' CP Position:',3f20.14)") CPpos(:,indcp) 1124 write(*,"(' CP type: ',a)") CPtyp2lab(CPtype(indcp)) 1125 call showptprop(CPpos(1,indcp),CPpos(2,indcp),CPpos(3,indcp),ifunctopo,6) 1126 else 1127 write(*,*) "Error: Invalid input" 1128 end if 1129 end if 1130!8888888888888888888 1131 else if (isel==8) then 1132 numpathold=numpath 1133 do i=1,numcp 1134 if (CPtype(i)==2) call findpath(i,1,ifunctopo) 1135 end do 1136 write(*,"(' Totally found',i6,' new paths')") numpath-numpathold 1137!9999999999999999999 1138 else if (isel==9) then 1139 numpathold=numpath 1140 do i=1,numcp 1141 if (CPtype(i)==3) call findpath(i,2,ifunctopo) 1142 end do 1143 write(*,"(' Totally found',i6,' new paths')") numpath-numpathold 1144!10 10 10 10 10 10 10 1145 else if (isel==10.and.count(CPtype(1:numcp)==2)==0) then 1146 write(*,*) "Error: You have to find at least one (3,-1) critical point" 1147 else if (isel==10) then 1148 do while(.true.) 1149 write(*,*) 1150 write(*,*) "Generate or delete interbasin surface for which (3,-1)?" 1151 write(*,*) "e.g. 5 means generate the surface from the (3,-1) with index of 5" 1152 if (numbassurf> 0) write(*,*) " -3 means delete the surface from the (3,-1) with index of 3" 1153 if (numbassurf> 0) write(*,*) "If input 0, surfaces from all (3,-1) will be deleted" 1154 if (numbassurf==0) write(*,*) "If input 0, surfaces from all (3,-1) will be generated" 1155 if (numbassurf> 0) write(*,*) "To list all generated surfaces, input the letter ""l""" 1156 if (numbassurf> 0) write(*,*) "To export the surfaces from the (3,-1) with index of 4, input ""o 4""" 1157 write(*,*) "To return, input ""q""" 1158 read(*,"(a)") c200 1159 1160 if (c200(1:1)=='q') then 1161 exit 1162 else if (c200(1:1)=='l') then 1163 do isurf=1,numbassurf 1164 do icp=1,numcp 1165 if (cp2surf(icp)==isurf) exit 1166 end do 1167 write(*,"('Index of surface:',i8,' Index of corresponding (3,-1):',i8)") isurf,icp 1168 end do 1169 else if (c200(1:1)=='o') then 1170 read(c200(3:),*) icp 1171 if (icp>numcp.or.icp<=0) then 1172 write(*,*) "The index of the surface is nonexisted" 1173 else if (cp2surf(icp)==0) then 1174 write(*,*) "This CP is not (3,-1), input again" 1175 else 1176 open(10,file="surpath.txt",status="replace") 1177 do ipath=1,nsurfpathpercp 1178 write(10,"('Path',i8)") ipath 1179 do ipt=1,nsurfpt 1180 write(10,"(i6,3f14.8)") ipt,bassurpath(:,ipt,ipath,cp2surf(icp)) 1181 end do 1182 end do 1183 close(10) 1184 write(*,"(a)") "The coordinates of the paths of the surface have been exported to surpath.txt in current folder" 1185 end if 1186 else 1187 read(c200,*) isel2 1188 if (isel2==0) then 1189 if (numbassurf>0) then 1190 numbassurf=0 1191 cp2surf=0 !If cps2surf(i)==0, means the (3,-1) with total index of i hasn't been given surface 1192 deallocate(bassurpath) 1193 else if (numbassurf==0) then !Generate all interbasin surfaces 1194 numbassurf=count(CPtype(1:numcp)==2) 1195 allocate(bassurpath(3,nsurfpt,nsurfpathpercp,numbassurf)) 1196 isurf=1 1197 write(*,"(i8,' surfaces will be generated, please wait...')") numbassurf 1198 do icp=1,numcp 1199 if (CPtype(icp)==2) then 1200 cp2surf(icp)=isurf 1201 call genbassurf(icp,isurf,ifunctopo) 1202 isurf=isurf+1 1203 write(*,"(a,i6)") " Finished the surface generation from (3,-1) with index of",icp 1204 end if 1205 end do 1206 end if 1207 else if (abs(isel2)>numcp) then 1208 write(*,*) "Error: This CP is nonexisted" 1209 else if (CPtype(abs(isel2))/=2) then 1210 write(*,*) "This CP is not (3,-1), input again" 1211 else if (isel2>0) then !Add a surface 1212 if (cp2surf(isel2)/=0) then 1213 write(*,*) "The interbasin surface has already been generated" 1214 else 1215 write(*,*) "Please wait..." 1216 if (numbassurf>0) then 1217 allocate(bassurpathtmp(3,nsurfpt,nsurfpathpercp,numbassurf)) 1218 bassurpathtmp=bassurpath 1219 deallocate(bassurpath) 1220 numbassurf=numbassurf+1 1221 allocate(bassurpath(3,nsurfpt,nsurfpathpercp,numbassurf)) 1222 bassurpath(:,:,:,1:numbassurf-1)=bassurpathtmp(:,:,:,:) 1223 deallocate(bassurpathtmp) 1224 else if (numbassurf==0) then 1225 numbassurf=numbassurf+1 1226 allocate(bassurpath(3,nsurfpt,nsurfpathpercp,numbassurf)) 1227 end if 1228 call genbassurf(isel2,numbassurf,ifunctopo) 1229 cp2surf(isel2)=numbassurf 1230 write(*,*) "Done!" 1231 end if 1232 else if (isel2<0) then !Delete a surface 1233 idelsurf=cp2surf(abs(isel2)) 1234 if (idelsurf==0) then 1235 write(*,*) "Surface has not been generated from this CP" 1236 else 1237 allocate(bassurpathtmp(3,nsurfpt,nsurfpathpercp,numbassurf)) 1238 bassurpathtmp=bassurpath 1239 deallocate(bassurpath) 1240 numbassurf=numbassurf-1 1241 allocate(bassurpath(3,nsurfpt,nsurfpathpercp,numbassurf)) 1242 bassurpath(:,:,:,1:idelsurf-1)=bassurpathtmp(:,:,:,1:idelsurf-1) 1243 bassurpath(:,:,:,idelsurf:)=bassurpathtmp(:,:,:,idelsurf+1:) 1244 deallocate(bassurpathtmp) 1245 cp2surf(abs(isel2))=0 1246 where(cp2surf>idelsurf) cp2surf=cp2surf-1 1247 end if 1248 end if 1249 end if 1250 end do 1251 else if (isel==20.and.ifunctopo==1) then 1252 do while(.true.) 1253 write(*,*) "Input the indices of the CPs in the ring, e.g. 22,23,25,28,32,11" 1254 write(*,*) "(Input q can exit)" 1255 read(*,"(a)") c200 1256 if (c200(1:1)=='q'.or.c200(1:1)=='Q') exit 1257 call str2arr(c200,nshanaromat,shanCPind) 1258 allocate(shanCPrho(nshanaromat)) 1259 totdens=0D0 1260 do ishan=1,nshanaromat 1261 shanCPrho(ishan)=fdens(CPpos(1,shanCPind(ishan)),CPpos(2,shanCPind(ishan)),CPpos(3,shanCPind(ishan))) 1262 end do 1263 totrho=sum(shanCPrho(1:nshanaromat)) 1264 shant=0D0 1265 do ishan=1,nshanaromat 1266 shanentropy=-shanCPrho(ishan)/totrho*log(shanCPrho(ishan)/totrho) 1267 write(*,"(' Electron density at CP',i3,':',f15.10,' Local entropy:',f15.10)") shanCPind(ishan),shanCPrho(ishan),shanentropy 1268 shant=shant+shanentropy 1269 end do 1270 shanmax=log(dfloat(nshanaromat)) 1271 write(*,"(' Total electron density:',f15.10)") totrho 1272 write(*,"(' Total Shannon entropy:',f15.10)") shant 1273 write(*,"(' Expected maximum Shannon entropy:',f15.10)") shanmax 1274 write(*,*) 1275 write(*,"(' Shannon aromaticity index:',f16.10)") shanmax-shant 1276 deallocate(shanCPrho) 1277 end do 1278 else if (isel==21.and.ifunctopo==1) then 1279 write(*,*) "Input the coordinate, e.g. 2.0,2.4,1.1 or input indices of a CP, e.g. 4" 1280 read(*,"(a)") c200 1281 if ( index(c200,',')==0 .and. index(trim(c200),' ')==0 ) then 1282 read(c200,*) ithisCP 1283 tmpx=CPpos(1,ithisCP) 1284 tmpy=CPpos(2,ithisCP) 1285 tmpz=CPpos(3,ithisCP) 1286 else 1287 read(c200,*) tmpx,tmpy,tmpz 1288 write(*,*) "The coordinate you inputted is in which unit? 1=Bohr 2=Angstrom" 1289 read(*,*) iunit 1290 if (iunit==2) then 1291 tmpx=tmpx/b2a 1292 tmpy=tmpy/b2a 1293 tmpz=tmpz/b2a 1294 end if 1295 end if 1296 write(*,*) "Input indices of three atoms to define a plane, e.g. 3,4,9" 1297 read(*,*) iatm1,iatm2,iatm3 1298 call pointABCD(a(iatm1)%x,a(iatm1)%y,a(iatm1)%z,a(iatm2)%x,a(iatm2)%y,a(iatm2)%z,a(iatm3)%x,a(iatm3)%y,a(iatm3)%z,xnor,ynor,znor,rnouse) !Normal vector is (xnor,ynor,znor) 1299 facnorm=sqrt(xnor**2+ynor**2+znor**2) 1300 xnor=xnor/facnorm !Normalize normal vector, then (xnor,ynor,znor) is the unit vector normal to the plane defined by iatm1,iatm2,iatm3 1301 ynor=ynor/facnorm 1302 znor=znor/facnorm 1303 if (allocated(b)) then 1304 call gencalchessmat(2,1,tmpx,tmpy,tmpz,densvalue,gradtmp,hesstmp) 1305 densgrad=xnor*gradtmp(1)+ynor*gradtmp(2)+znor*gradtmp(3) 1306 denscurvature=xnor*xnor*hesstmp(1,1)+xnor*ynor*hesstmp(1,2)+xnor*znor*hesstmp(1,3)+& 1307 ynor*xnor*hesstmp(2,1)+ynor*ynor*hesstmp(2,2)+ynor*znor*hesstmp(2,3)+& 1308 znor*xnor*hesstmp(3,1)+znor*ynor*hesstmp(3,2)+znor*znor*hesstmp(3,3) 1309 write(*,"(' The unit normal vector is',3f14.8)") xnor,ynor,znor 1310 write(*,"(' Electron density is ',f30.10)") densvalue 1311 write(*,"(' Electron density gradient is ',f30.10)") densgrad 1312 write(*,"(' Electron density curvature is',f30.10)") denscurvature 1313 end if 1314 write(*,*) 1315 write(*,"(a)") " BTW: The X,Y,Z coordinate (row) of current point, the points below and above 1 Angstrom of the plane from current point, respectively (in Angstrom)." 1316 write(*,"(3f16.10)") tmpx*b2a,tmpy*b2a,tmpz*b2a 1317 write(*,"(3f16.10)") (tmpx-xnor/b2a)*b2a,(tmpy-ynor/b2a)*b2a,(tmpz-znor/b2a)*b2a 1318 write(*,"(3f16.10)") (tmpx+xnor/b2a)*b2a,(tmpy+ynor/b2a)*b2a,(tmpz+znor/b2a)*b2a 1319 write(*,*) 1320 end if 1321end do 1322end subroutine 1323 1324 1325!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1326!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1327 1328!!!--------- Generate interbasin surface from (3,-1) 1329subroutine genbassurf(ithisCP,ithissurf,ifunc) 1330use defvar 1331use function 1332use util 1333use topo 1334implicit real*8(a-h,o-z) 1335integer ifunc,ithisCP,ithissurf 1336real*8 initstpsize,hess(3,3),grad(3),eigvecmat(3,3),eigval(3),basvec1(3),basvec2(3),k1(3),k2(3) 1337initstpsize=surfpathstpsiz/4D0 !smaller than pathstepsize(0.02) 1338 1339call gencalchessmat(2,ifunc,CPpos(1,ithisCP),CPpos(2,ithisCP),CPpos(3,ithisCP),value,grad,hess) 1340call diagmat(hess,eigvecmat,eigval,300,1D-12) 1341!Generate normalized basis vectors from two negative eignvectors of hessian matrix 1342itmp=0 1343do i=1,3 1344 if (eigval(i)<0) then 1345 if (itmp==0) then 1346 rnorm=dsqrt(eigvecmat(1,i)**2+eigvecmat(2,i)**2+eigvecmat(3,i)**2) 1347 basvec1=eigvecmat(:,i)/rnorm 1348 itmp=1 1349 else if (itmp==1) then 1350 rnorm=dsqrt(eigvecmat(1,i)**2+eigvecmat(2,i)**2+eigvecmat(3,i)**2) 1351 basvec2=eigvecmat(:,i)/rnorm 1352 end if 1353 end if 1354end do 1355!Use the two basis vectors to generate a circle of initial points around (3,-1) 1356angstp=360D0/nsurfpathpercp 1357do i=1,nsurfpathpercp 1358 ang=(i-1)*angstp/180D0*pi !Convert to arc unit 1359 bassurpath(:,1,i,ithissurf)=initstpsize*( sin(ang)*basvec1(:)+cos(ang)*basvec2(:) ) + CPpos(:,ithisCP) 1360end do 1361 1362!Generate gradient path from each initial point to comprise interbasin surface 1363do ipath=1,nsurfpathpercp 1364 do ipt=2,nsurfpt 1365 !Move point, RK2 method. Only calculate function value and gradient 1366 xtmp=bassurpath(1,ipt-1,ipath,ithissurf) 1367 ytmp=bassurpath(2,ipt-1,ipath,ithissurf) 1368 ztmp=bassurpath(3,ipt-1,ipath,ithissurf) 1369 call gencalchessmat(1,ifunc,xtmp,ytmp,ztmp,value,grad,hess) 1370 k1=grad/dsqrt(sum(grad**2)) 1371 call gencalchessmat(1,ifunc,xtmp+surfpathstpsiz/2*k1(1),ytmp+surfpathstpsiz/2*k1(2),ztmp+surfpathstpsiz/2*k1(3),value,grad,hess) 1372 k2=grad/dsqrt(sum(grad**2)) 1373 bassurpath(:,ipt,ipath,ithissurf)=bassurpath(:,ipt-1,ipath,ithissurf)-surfpathstpsiz*k2 1374 end do 1375end do 1376end subroutine 1377 1378 1379 1380!!--------- Generate interbasin path from (3,-1) on a given plane 1381!ifunc=which real space function, ithisCP: The total index of (3,-1), ithispath: Store to which slot of ple3n1path array 1382subroutine gen3n1plepath(ifunc,ithisCP,ithispath) 1383use defvar 1384use topo 1385use function 1386use util 1387implicit real*8(a-h,o-z) 1388integer ifunc,ithisCP,ithispath 1389real*8 initstpsize,hess(3,3),grad(3),eigvecmat(3,3),eigval(3),initvec(3),plenormvec(3),k1(3),k2(3) 1390initstpsize=ple3n1pathstpsiz/4D0 !smaller than pathstepsize(0.02) 1391 1392call gencalchessmat(2,ifunc,CPpos(1,ithisCP),CPpos(2,ithisCP),CPpos(3,ithisCP),value,grad,hess) 1393call diagmat(hess,eigvecmat,eigval,300,1D-12) 1394do iposvec=1,3 1395 if (eigval(iposvec)>0) exit 1396end do 1397if (plesel==1) then !XY 1398 plenormvec(1:2)=0D0 1399 plenormvec(3)=1D0 1400else if (plesel==2) then !XZ 1401 plenormvec(1:3)=0D0 1402 plenormvec(2)=1D0 1403else if (plesel==3) then !YZ 1404 plenormvec(2:3)=0D0 1405 plenormvec(1)=1D0 1406else 1407 call pointABCD(a1x,a1y,a1z,a2x,a2y,a2z,a3x,a3y,a3z,plenormvec(1),plenormvec(2),plenormvec(3),tmp) 1408end if 1409 1410call vecprod(plenormvec(1),plenormvec(2),plenormvec(3), eigvecmat(1,iposvec),eigvecmat(2,iposvec),eigvecmat(3,iposvec), initvec(1),initvec(2),initvec(3)) 1411rnorm=dsqrt(sum(initvec**2)) 1412initvec=initvec/rnorm 1413ple3n1path(:,1,1,ithispath)=CPpos(:,ithisCP)+initvec(:)*initstpsize 1414ple3n1path(:,1,2,ithispath)=CPpos(:,ithisCP)-initvec(:)*initstpsize 1415 1416!Generate gradient path from the given (3,-1) to comprise interbasin path on the plane 1417do idir=1,2 1418 do ipt=2,n3n1plept 1419 !Move point, RK2 method. Only calculate function value and gradient, don't calculate Hessian for saving time 1420 xtmp=ple3n1path(1,ipt-1,idir,ithispath) 1421 ytmp=ple3n1path(2,ipt-1,idir,ithispath) 1422 ztmp=ple3n1path(3,ipt-1,idir,ithispath) 1423 call gencalchessmat(1,ifunc,xtmp,ytmp,ztmp,value,grad,hess) 1424 k1=grad/dsqrt(sum(grad**2)) 1425 call gencalchessmat(1,ifunc,xtmp+surfpathstpsiz/2*k1(1),ytmp+surfpathstpsiz/2*k1(2),ztmp+surfpathstpsiz/2*k1(3),value,grad,hess) 1426 k2=grad/dsqrt(sum(grad**2)) 1427 ple3n1path(:,ipt,idir,ithispath)=ple3n1path(:,ipt-1,idir,ithispath)-surfpathstpsiz*k2 1428 end do 1429end do 1430end subroutine 1431 1432 1433 1434!!!--------- Find path from a critical point at x,y,z 1435!itype=1: from (3,-1) to (3,-3) =2: from (3,+1) to (3,+3) =3: between (3,+1) and (3,-1) 1436subroutine findpath(ithisCP,itype,ifunc) 1437use topo 1438use function 1439use util 1440implicit real*8(a-h,o-z) 1441integer ifunc,ithisCP,foundind(npathtry) !foundind records that in pathtmp, which saved newly found path, =1 means saved 1442real*8 grad(3),hess(3,3),k1(3),k2(3) 1443real*8 eigvecmat(3,3),eigval(3),pathtmp(3,maxpathpt,npathtry) !Will maximally try to find npathtry paths 1444real*8,allocatable :: tmparr(:,:,:) 1445!Note: Each time invoke this routine, for itype=1 and 2, search two times; for itype=3, search npathtry times 1446!The new path first store to pathtmp, then enlarge topopath and pass the new path information to it 1447 1448foundind=0 1449noldpath=numpath 1450ipath=0 1451!Determine eigenvalue and eigenvector of Hessian at initial point 1452call gencalchessmat(2,ifunc,CPpos(1,ithisCP),CPpos(2,ithisCP),CPpos(3,ithisCP),value,grad,hess) 1453call diagmat(hess,eigvecmat,eigval,300,1D-12) 1454 1455if (itype==1.or.itype==2) then 1456 pathtmp(:,1,1)=CPpos(:,ithisCP) !Set first point as input coordinate 1457 pathtmp(:,1,2)=CPpos(:,ithisCP) 1458 if (itype==1) then 1459 do iposi=1,3 1460 if (eigval(iposi)>0) exit !Find positive eigenvalue of Hessian of (3,-1) 1461 end do 1462 else 1463 do iposi=1,3 1464 if (eigval(iposi)<0) exit !Find negative eigenvalue of Hessian of (3,+1) 1465 end do 1466 end if 1467iterdir: do idir=1,2 1468 if (idir==1) write(*,"(' Go forward from CP: ',i6,1x,a,' Position:',3f12.6)") ithisCP,CPtyp2lab(CPtype(ithisCP)),CPpos(1:3,ithisCP) 1469 if (idir==2) write(*,"(' Go backward from CP:',i6,1x,a,' Position:',3f12.6)") ithisCP,CPtyp2lab(CPtype(ithisCP)),CPpos(1:3,ithisCP) 1470 posvecnorm=dsqrt(sum(eigvecmat(:,iposi)**2)) 1471 if (idir==1) pathtmp(:,2,idir)=pathtmp(:,1,idir)+pathstepsize*eigvecmat(:,iposi)/posvecnorm !Move forwards along eigenvector with positive eigenvalue 1472 if (idir==2) pathtmp(:,2,idir)=pathtmp(:,1,idir)-pathstepsize*eigvecmat(:,iposi)/posvecnorm !Move backwards along eigenvector with positive eigenvalue 1473 !Check if the path has already presented by comparing corresponding first two points 1474 if (allocated(topopath)) then 1475 do ickpath=1,noldpath 1476! write(*,*) ickpath,numpath,size(topopath,3) 1477 if ( dsqrt(sum( (topopath(:,1,ickpath)-pathtmp(:,1,idir))**2 ))<0.01D0 & 1478 .and.dsqrt(sum( (topopath(:,2,ickpath)-pathtmp(:,2,idir))**2 ))<0.01D0) then 1479 write(*,*) "The path has already presented, skip search" 1480 write(*,*) 1481 cycle iterdir 1482 end if 1483 end do 1484 end if 1485 1486iterpt: do ipt=2,maxpathpt 1487 !Check if the distance between current point and any other CP is smaller than threshold (discritpathfin) 1488 do icp=1,numcp 1489 if (icp==ithisCP) cycle 1490 distcp=dsqrt(sum( (pathtmp(:,ipt,idir)-CPpos(:,icp))**2 )) 1491 if (distcp<discritpathfin) then 1492 numpath=numpath+1 1493 foundind(idir)=1 1494 pathnumpt(numpath)=ipt 1495 write(*,"(a,i6,a,f8.4)") " Found new path after",ipt," iterations, path length:",(ipt-1)*pathstepsize 1496 write(*,"(' Reached CP',i6,1x,a,' Position:',3f12.6)") icp,CPtyp2lab(CPtype(icp)),CPpos(:,icp) 1497! do i=1,ipt 1498! write(*,"(i6,3f12.6)") i,pathtmp(:,i,idir) 1499! end do 1500 exit iterpt 1501 end if 1502 end do 1503 if (ipt==maxpathpt) then 1504 write(*,*) "Search failed, exceeded upper limit of number of iterations" 1505 exit 1506 end if 1507 1508 !Move point, RK2 method. Only calculate function value and gradient 1509 xtmp=pathtmp(1,ipt,idir) 1510 ytmp=pathtmp(2,ipt,idir) 1511 ztmp=pathtmp(3,ipt,idir) 1512 call gencalchessmat(1,ifunc,xtmp,ytmp,ztmp,value,grad,hess) 1513 k1=grad/dsqrt(sum(grad**2)) 1514 call gencalchessmat(1,ifunc,xtmp+pathstepsize/2*k1(1),ytmp+pathstepsize/2*k1(2),ztmp+pathstepsize/2*k1(3),value,grad,hess) 1515 k2=grad/dsqrt(sum(grad**2)) 1516 if (itype==1) pathtmp(:,ipt+1,idir)=pathtmp(:,ipt,idir)+pathstepsize*k2 1517 if (itype==2) pathtmp(:,ipt+1,idir)=pathtmp(:,ipt,idir)-pathstepsize*k2 1518 end do iterpt 1519 1520 write(*,*) 1521 end do iterdir 1522else if (itype==3) then 1523end if 1524 1525!Enlarge old topopath and then add the path found at this time to it 1526numnewpath=numpath-noldpath !How many new path found in this time 1527if (numnewpath/=0) then 1528 itmp=1 1529 if (allocated(topopath)) then 1530 allocate(tmparr(3,maxpathpt,noldpath)) 1531 tmparr=topopath 1532 deallocate(topopath) 1533 allocate(topopath(3,maxpathpt,numpath)) 1534 topopath(:,:,1:noldpath)=tmparr 1535 do i=1,npathtry !The number of foundind(i)==1 equals numnewpath 1536 if (foundind(i)==1) then 1537 topopath(:,:,noldpath+itmp)=pathtmp(:,:,i) 1538 itmp=itmp+1 1539 end if 1540 end do 1541 deallocate(tmparr) 1542 else !This first time invoke this routine 1543 allocate(topopath(3,maxpathpt,numnewpath)) 1544 do i=1,npathtry 1545 if (foundind(i)==1) then 1546 topopath(:,:,itmp)=pathtmp(:,:,i) 1547 itmp=itmp+1 1548 end if 1549 end do 1550 end if 1551end if 1552end subroutine 1553 1554 1555 1556!!--------- Determine path type and connected which two CPs, 0 means unknown (not connected a found CP) 1557!ipathtype =0: other =1: (3,-1)->(3,-3) =2: (3,+1)->(3,+3) =3: (3,-1)<-->(3,+1) 1558subroutine path_cp(ipath,icp1,icp2,ipathtype) 1559use topo 1560implicit real*8(a-h,o-z) 1561integer ipath,icp1,icp2 1562iunknown=0 1563do icp1=1,numcp 1564 if (sum( (topopath(:,1,ipath)-CPpos(:,icp1))**2 )<discritpathfin**2) exit !Test the first point in the path 1565 if (icp1==numcp) iunknown=1 1566end do 1567if (iunknown==1) icp1=0 1568iunknown=0 1569do icp2=1,numcp 1570 if (sum( (topopath(:,pathnumpt(ipath),ipath)-CPpos(:,icp2))**2 )<discritpathfin**2) exit !Test the last point 1571 if (icp2==numcp) iunknown=0 1572end do 1573if (iunknown==1) icp2=0 1574 1575if (icp1==0.or.icp2==0) then 1576 ipathtype=0 1577else 1578 icp1type=CPtype(icp1) 1579 icp2type=CPtype(icp2) 1580 if (icp1type==2.and.icp2type==1) then 1581 ipathtype=1 1582 else if (icp1type==3.and.icp2type==4) then 1583 ipathtype=2 1584 else if ( (icp1type==2.and.icp2type==3).or.(icp1type==3.and.icp2type==2) ) then 1585 ipathtype=3 1586 else 1587 ipathtype=0 1588 end if 1589end if 1590end subroutine 1591 1592 1593 1594!!!----------- Find critical points from initial guess at X,Y,Z 1595!ifunc is index of real space functions 1596!If ilowcrit==1, use lower critiera, because for heavy atoms, cusp of electron density at nuclear position is very sharp hence hard to locate by default criteria 1597!ishowsearchlevel=0/1/2: Print none/some detail/all detail. Notice that in parallel mode, the outputted details are messed up 1598subroutine findcp(x,y,z,ifunc,ilowcrit) 1599use topo 1600use function 1601use util 1602implicit real*8(a-h,o-z) 1603integer ifunc 1604integer ilowcrit 1605real*8 x,y,z 1606real*8 coord(3,1),grad(3,1),hess(3,3),disp(3,1) 1607real*8 eigvecmat(3,3),eigval(3) !,tmpmat(3,3) 1608if (ilowcrit==1) then 1609 realdispconv=dispconv*10000D0 1610 realgradconv=gradconv*10000D0 1611else !For most cases use default criteria 1612 realdispconv=dispconv 1613 realgradconv=gradconv 1614end if 1615coord(1,1)=x 1616coord(2,1)=y 1617coord(3,1)=z 1618if (ishowsearchlevel>=1) write(*,"(' Starting point:',3f12.6)") coord(1:3,1) 1619 1620do i=1,topomaxcyc 1621 call gencalchessmat(2,ifunc,coord(1,1),coord(2,1),coord(3,1),value,grad(1:3,1),hess) 1622! call showmatgau(hess) 1623! write(*,*) detmat(hess) 1624 singulartest=abs(detmat(hess)) 1625 if (singulartest<singularcrit) then 1626 if (ishowsearchlevel>=1) then 1627 write(*,*) "Hessian matrix is singular at current position, stop iteration" 1628 write(*,"(' Absolute of determinant of Hessian matrix:',E18.10)") singulartest 1629 write(*,"(' Criterion for detecting singular:',E18.10)") singularcrit 1630 end if 1631 exit 1632 end if 1633 disp=-matmul(invmat(hess,3),grad) 1634 coord=coord+CPstepscale*disp 1635 disperr=dsqrt(sum(disp**2)) 1636 graderr=dsqrt(sum(grad**2)) 1637 1638 if (ishowsearchlevel==2) then 1639 write(*,"(/,' Step',i5,' Function Value:',f18.10)") i,value 1640 write(*,"(' Displacement vector:',3f18.10)") disp 1641 write(*,"(' Current coordinate :',3f18.10)") coord 1642 write(*,"(' Current gradient :',3E18.10)") grad 1643 write(*,"(' Norm of displacement:',E18.8,' Norm of gradient:',E18.8)") disperr,graderr 1644 write(*,"(' Goal: |disp|<',E18.8,' |Grad|<',E18.8)") realdispconv,realgradconv 1645 end if 1646! tmpmat=hess 1647! call diagmat(tmpmat,eigvecmat,eigval,300,1D-12) 1648! write(*,"('Eigenvalue of Hessian : ',3D16.8)") eigval 1649 1650 if (disperr<realdispconv.and.graderr<realgradconv) then 1651 if (ishowsearchlevel>=1) write(*,"(' After',i6,' iterations')") i 1652 if (ishowsearchlevel==2) write(*,*) "====================== Iteration ended ======================" 1653 inewcp=1 1654 do icp=1,numcp 1655 r=dsqrt(sum( (coord(:,1)-CPpos(:,icp))**2 )) 1656 if (r<=minicpdis) then 1657 if (ishowsearchlevel>=1) write(*,"(a,i6,a)") " This CP is too close to CP",icp,", ignored..." 1658 inewcp=0 1659 exit 1660 end if 1661 end do 1662 if (CPsearchlow/=CPsearchhigh.and.(value<CPsearchlow.or.value>CPsearchhigh)) then 1663 write(*,"(a,1PD12.5,a)") " The value of this CP is ",value,", which exceeded user-defined range and thus ignored" 1664 inewcp=0 1665 end if 1666 if (inewcp==1) then 1667!$OMP CRITICAL 1668 numcp=numcp+1 1669 CPpos(:,numcp)=coord(:,1) 1670 call diagmat(hess,eigvecmat,eigval,300,1D-15) 1671! call diagsymat(hess,eigvecmat,eigval,idiagok) 1672! if (idiagok/=0) write(*,*) "Diagonization of Hessian matrix failed" 1673 igt0=count(eigval>0) 1674 if (igt0==3) then 1675 if (ishowsearchlevel>=1) write(*,"(' Found new (3,+3) at',3f15.10)") coord 1676 CPtype(numcp)=4 1677 else if (igt0==2) then 1678 if (ishowsearchlevel>=1) write(*,"(' Found new (3,+1) at',3f15.10)") coord 1679 CPtype(numcp)=3 1680 else if (igt0==1) then 1681 if (ishowsearchlevel>=1) write(*,"(' Found new (3,-1) at',3f15.10)") coord 1682 CPtype(numcp)=2 1683 call sort(eigval) 1684 if (ishowsearchlevel>=1) write(*,"(' Bond ellipticity is',f15.10)") eigval(1)/eigval(2)-1.0D0 1685 else if (igt0==0) then 1686 if (ishowsearchlevel>=1) write(*,"(' Found new (3,-3) at',3f15.10)") coord 1687 CPtype(numcp)=1 1688 end if 1689 if (ishowsearchlevel>=1) write(*,"(' Eigenvalue:',3f20.10)") eigval 1690!$OMP end CRITICAL 1691 end if 1692 exit 1693 end if 1694 if (i==topomaxcyc.and.(ishowsearchlevel>=1)) write(*,*) "!! Exceeded maximum cycle until find stationary point !!" 1695end do 1696if (ishowsearchlevel>=1) write(*,*) 1697end subroutine 1698 1699 1700!Sort newly found CPs according to coordinates 1701!numcpoldp1: The number of CPs before this search + 1 1702subroutine sortCP(numcpoldp1) 1703use topo 1704implicit real*8(a-h,o-z) 1705integer numcpoldp1,typetmp 1706real*8 tmparr(3) 1707do itmp=numcpoldp1,numcp 1708 do jtmp=itmp+1,numcp 1709 tmpvali=CPpos(1,itmp)*0.234134D0+CPpos(2,itmp)*1.9837322D0-CPpos(3,itmp)*0.5413578924D0 !Use three random number to generate unique code 1710 tmpvalj=CPpos(1,jtmp)*0.234134D0+CPpos(2,jtmp)*1.9837322D0-CPpos(3,jtmp)*0.5413578924D0 1711 if (tmpvali>tmpvalj) then 1712 typetmp=CPtype(itmp) 1713 CPtype(itmp)=CPtype(jtmp) 1714 CPtype(jtmp)=typetmp 1715 tmparr=CPpos(:,itmp) 1716 CPpos(:,itmp)=CPpos(:,jtmp) 1717 CPpos(:,jtmp)=tmparr 1718 end if 1719 end do 1720end do 1721end subroutine 1722