1subroutine bondana 2use defvar 3use util 4implicit real*8 (a-h,o-z) 5character c2000tmp*2000 6do while(.true.) 7 write(*,*) " ================ Bond order analysis ===============" 8 if (allocated(CObasa)) then 9 if (allocated(frag1)) then 10 write(*,*) "-1 Redefine fragment 1 and 2 for option 1,3,4,7,8" 11 else 12 write(*,*) "-1 Define fragment 1 and 2 for option 1,3,4,7,8 (to be defined)" 13 end if 14 end if 15 write(*,*) "0 Return" 16 write(*,*) "1 Mayer bond order analysis" 17 write(*,*) "2 Multicenter bond order analysis (up to 12 centers)" 18 write(*,*) "-2 Multicenter bond order analysis in NAO basis (up to 12 centers)" 19! write(*,*) "-3 Multicenter bond order analysis in Lowdin orthogonalized basis (up to 12 centers)" !Can be used, but not very meaningful, so not shown 20 write(*,*) "3 Wiberg bond order analysis in Lowdin orthogonalized basis" 21 write(*,*) "4 Mulliken bond order analysis" 22 write(*,*) "5 Decompose Mulliken bond order between two atoms to orbital contributions" 23 write(*,*) "6 Orbital occupancy-perturbed Mayer bond order" 24 write(*,*) "7 Fuzzy bond order analysis" 25 write(*,*) "8 Laplacian bond order" 26 write(*,*) "9 Decompose Wiberg bond order in NAO basis as atomic orbital pair contribution" 27 read(*,*) ibondana 28 if (.not.allocated(CObasa).and.(ibondana>=1.and.ibondana<=6)) then 29 write(*,"(a)") " ERROR: The input file you used does not contain basis function information! Please check Section 2.5 of the manual for explanation" 30 return 31 else if (.not.allocated(b).and.(ibondana==7.or.ibondana==8)) then 32 write(*,"(a)") " ERROR: The input file you used does not contain GTF information! Please check Section 2.5 of the manual for explanation" 33 return 34 end if 35 36 if (ibondana==-1) then 37 !Define frag1, the size just accomodates content 38 if (allocated(frag1)) then 39 write(*,*) "Atoms in current fragment 1:" 40 write(*,"(13i6)") frag1 41 write(*,"(a)") " Input 0 to keep unchanged, or redefine fragment 1, e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment 1" 42 else 43 write(*,"(a)") " Input atomic indices to define fragment 1. e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment 1" 44 end if 45 read(*,"(a)") c2000tmp 46 if (c2000tmp(1:1)/='0') then 47 if (allocated(frag1)) deallocate(frag1) 48 call str2arr(c2000tmp,nfragatm) 49 allocate(frag1(nfragatm)) 50 call str2arr(c2000tmp,nfragatm,frag1) 51 end if 52 !Define frag2, the size just accomodates content 53 if (allocated(frag2)) then 54 write(*,*) "Atoms in current fragment 2:" 55 write(*,"(13i6)") frag2 56 write(*,"(a)") " Input 0 to keep unchanged, or redefine fragment 2, e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute fragment 2" 57 else 58 write(*,"(a)") " Input atomic indices to define fragment 2. e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment 2" 59 end if 60 read(*,"(a)") c2000tmp 61 if (c2000tmp(1:1)/='0') then 62 if (allocated(frag2)) deallocate(frag2) 63 call str2arr(c2000tmp,nfragatm) 64 allocate(frag2(nfragatm)) 65 call str2arr(c2000tmp,nfragatm,frag2) 66 end if 67 if (any(frag1>ncenter).or.any(frag2>ncenter)) then 68 write(*,*) "Error: Some atomic indices exceeded valid range! Please define again" 69 write(*,*) 70 deallocate(frag1,frag2) 71 cycle 72 end if 73 write(*,*) "Setting is saved" 74 write(*,*) "Now the atoms in fragment 1 are" 75 write(*,"(13i6)") frag1 76 write(*,*) "Now the atoms in fragment 2 are" 77 write(*,"(13i6)") frag2 78 if (any(frag1<=0).or.any(frag1>ncenter).or.any(frag2<=0).or.any(frag2>ncenter)) write(*,*) "Warning: Indices of some atoms exceed valid range! Please redefine fragment" 79 do i=1,size(frag1) 80 if (any(frag2==frag1(i))) then 81 write(*,"(a)") "Warning: Indices of some atoms are duplicated in the two fragments! Please redefine them" 82 exit 83 end if 84 end do 85 write(*,*) 86 87 else if (ibondana==0) then 88 if (allocated(frag1)) deallocate(frag1) 89 if (allocated(frag2)) deallocate(frag2) 90 exit 91 else if (ibondana==1) then 92 write(*,*) "Please wait..." 93 call mayerbndord 94 else if (ibondana==2) then 95 call multicenter 96 else if (ibondana==-2) then 97 call bndordNAO(1) 98 else if (ibondana==3.or.ibondana==-3) then 99 !In symmortho the density matrix, cobasa/b and Sbas will change, so backup them 100 if (allocated(Cobasb)) then !Open-shell 101 allocate(Cobasa_org(nbasis,nmo/2),Cobasb_org(nbasis,nmo/2),Sbas_org(nbasis,nbasis)) 102 Cobasb_org=Cobasb 103 else 104 allocate(Cobasa_org(nbasis,nmo),Sbas_org(nbasis,nbasis)) 105 end if 106 Cobasa_org=Cobasa 107 Sbas_org=Sbas 108 write(*,*) "Performing Lowdin orthogonalization..." 109 call symmortho 110 if (ibondana==3) then 111 write(*,*) "Calculating Wiberg bond order..." 112 call mayerbndord 113 else if (ibondana==-3) then 114 call multicenter 115 end if 116 write(*,*) "Regenerating original density matrix..." 117 write(*,*) 118 Cobasa=Cobasa_org 119 Sbas=Sbas_org 120 deallocate(Cobasa_org,Sbas_org) 121 if (allocated(Cobasb_org)) then 122 Cobasb=Cobasb_org 123 deallocate(Cobasb_org) 124 end if 125 call genP 126 else if (ibondana==4) then 127 call mullikenbndord 128 else if (ibondana==5) then 129 call decompMBO 130 else if (ibondana==6) then 131 call OrbPertMayer 132 else if (ibondana==7) then 133 call intatomspace(1) 134 else if (ibondana==8) then 135 call intatomspace(2) 136 else if (ibondana==9) then 137 call bndordNAO(2) 138 end if 139end do 140end subroutine 141 142 143!! ----------------- Mayer/Generalized Wiberg 2-c bond order analysis 144! Mayer bond order analysis and Generalized Wiberg bond order (GWBO) analysis 145! If Lowdin orthogonalization has been performed, that is carry out Wiberg bond order analysis in Lowdin orthogonalized basis 146! Note: For close-shell, two methods give the same result. For open-shell, Mayer bond order for all electrons is the sum of 147! alpha and beta bond order, while GWBO directly use total density matrix to generate 148! total bond order, the "Mayer bond order" in gaussian is actually GWBO! 149subroutine mayerbndord 150use defvar 151use util 152implicit real*8 (a-h,o-z) 153real*8 :: bndmata(ncenter,ncenter),bndmatb(ncenter,ncenter),bndmattot(ncenter,ncenter),& 154PSmata(nbasis,nbasis),PSmatb(nbasis,nbasis),PSmattot(nbasis,nbasis) 155character selectyn 156 157bndmata=0D0 158bndmatb=0D0 159bndmattot=0D0 160!Calculate total bond order for restricted close-shell wavefunction (for open-shell do GWBO, P=Palpha+Pbeta) 161PSmattot=matmul(Ptot,Sbas) 162do i=1,ncenter 163 do j=i+1,ncenter 164 accum=0D0 165 do ii=basstart(i),basend(i) 166 do jj=basstart(j),basend(j) 167 accum=accum+PSmattot(ii,jj)*PSmattot(jj,ii) 168 end do 169 end do 170 bndmattot(i,j)=accum 171 end do 172end do 173bndmattot=bndmattot+transpose(bndmattot) !Because we only filled one triangular region, copy it to another 174do i=1,ncenter 175 bndmattot(i,i)=sum(bndmattot(i,:)) 176end do 177 178if (wfntype==1.or.wfntype==2.or.wfntype==4) then 179 PSmata=matmul(Palpha,Sbas) 180 PSmatb=matmul(Pbeta,Sbas) 181 do i=1,ncenter 182 do j=i+1,ncenter 183 accuma=0D0 184 accumb=0D0 185 do ii=basstart(i),basend(i) 186 do jj=basstart(j),basend(j) 187 accuma=accuma+PSmata(ii,jj)*PSmata(jj,ii) 188 accumb=accumb+PSmatb(ii,jj)*PSmatb(jj,ii) 189 end do 190 end do 191 bndmata(i,j)=accuma 192 bndmatb(i,j)=accumb 193 end do 194 end do 195 bndmata=2*(bndmata+transpose(bndmata)) 196 bndmatb=2*(bndmatb+transpose(bndmatb)) 197 do i=1,ncenter 198 bndmata(i,i)=sum(bndmata(i,:)) 199 bndmatb(i,i)=sum(bndmatb(i,:)) 200 end do 201end if 202 203write(*,"(' The total bond order >=',f10.6)") bndordthres 204itmp=0 205if (wfntype==1.or.wfntype==2.or.wfntype==4) then 206 do i=1,ncenter 207 do j=i+1,ncenter 208 if (bndmata(i,j)+bndmatb(i,j)>=bndordthres) then 209 itmp=itmp+1 210 write(*,"(' #',i5,':',i5,a,i5,a,' Alpha: ',f10.6,' Beta:',f10.6,' Total:',f10.6)") & 211 itmp,i,'('//a(i)%name//')',j,'('//a(j)%name//')',bndmata(i,j),bndmatb(i,j),bndmata(i,j)+bndmatb(i,j) 212 end if 213 end do 214 end do 215 write(*,*) 216 write(*,"(' Bond order from mixed alpha&beta density matrix >=',f10.6)") bndordthres 217end if 218itmp=0 219do i=1,ncenter 220 do j=i+1,ncenter 221 if (bndmattot(i,j)>=bndordthres) then 222 itmp=itmp+1 223 write(*,"(' #',i5,':',5x,i5,a,i5,a,f14.8)") itmp,i,'('//a(i)%name//')',j,'('//a(j)%name//')',bndmattot(i,j) 224 end if 225 end do 226end do 227 228write(*,*) 229write(*,*) "Total valences and free valences defined by Mayer:" 230do i=1,ncenter 231 accum=0D0 232 accum2=0D0 233 do ii=basstart(i),basend(i) 234 accum=accum+2*PSmattot(ii,ii) 235 do jj=basstart(i),basend(i) 236 accum2=accum2+PSmattot(ii,jj)*PSmattot(jj,ii) 237 end do 238 end do 239 freeval=accum-accum2-(bndmata(i,i)+bndmatb(i,i)) 240 if (wfntype==0.or.wfntype==3) freeval=0D0 241 write(*,"(' Atom',i6,'(',a,') :',2f14.8)") i,a(i)%name,accum-accum2,freeval 242end do 243 244!Between fragment 245if (allocated(frag1)) then 246 bndordfraga=0 247 bndordfragb=0 248 bndordfragtot=0 249 do i=1,size(frag1) 250 do j=1,size(frag2) 251 bndordfraga=bndordfraga+bndmata(frag1(i),frag2(j)) 252 bndordfragb=bndordfragb+bndmatb(frag1(i),frag2(j)) 253 bndordfragtot=bndordfragtot+bndmattot(frag1(i),frag2(j)) 254 end do 255 end do 256 write(*,*) 257 if (wfntype==1.or.wfntype==2.or.wfntype==4) then 258 write(*,"(' The bond order between fragment 1 and 2:')") 259 write(*,"(' Alpha:',f10.6,' Beta:',f10.6,' Total:',f10.6,' Mixed Alpha&Beta:',f10.6)") bndordfraga,bndordfragb,bndordfraga+bndordfragb,bndordfragtot 260 else 261 write(*,"(' The bond order between fragment 1 and 2:',f12.6)") bndordfragtot 262 end if 263end if 264write(*,*) 265 266write(*,*) "If output bond order matrix to bndmat.txt in current folder? (y/n)" 267read(*,*) selectyn 268if (selectyn=='y'.or.selectyn=='Y') then 269 open(10,file="bndmat.txt",status="replace") 270 write(10,*) "Note: The diagonal elements are the sum of corresponding row elements" 271 if (wfntype==0.or.wfntype==3) then 272 call showmatgau(bndmattot,"Bond order matrix",0,"f14.8",10) 273 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 274 call showmatgau(bndmata,"Bond order matrix for alpha electrons",0,"f14.8",10) 275 call showmatgau(bndmatb,"Bond order matrix for beta electrons",0,"f14.8",10) 276 call showmatgau(bndmata+bndmatb,"Bond order matrix for all electrons",0,"f14.8",10) 277 call showmatgau(bndmattot,"Bond order matrix from mixed density",0,"f14.8",10) 278 end if 279 close(10) 280 write(*,*) "Result have been outputted to bndmat.txt in current folder" 281 write(*,*) 282end if 283end subroutine 284 285 286 287!! ----------- Multi-center bond order analysis 288subroutine multicenter 289use defvar 290use util 291implicit real*8 (a-h,o-z) 292real*8 :: PSmat(nbasis,nbasis),PSmata(nbasis,nbasis),PSmatb(nbasis,nbasis),PSmattot(nbasis,nbasis),maxbndord 293integer cenind(12),maxcenind(12) !maxcenind is used to store the combination that has the maximum bond order in automatical search 294integer itype 295character c1000tmp*1000 296 297write(*,*) "Please wait..." 298ntime=1 !Closed-shell 299PSmattot=matmul(Ptot,Sbas) 300if (wfntype==1.or.wfntype==2.or.wfntype==4) then !Open-shell 301 ntime=3 302 PSmata=matmul(Palpha,Sbas) 303 PSmatb=matmul(Pbeta,Sbas) 304end if 305 306do while(.true.) 307 write(*,*) 308 write(*,*) "Input atom indices, e.g. 3,4,7,8,10 (2~12 centers)" 309 write(*,*) "Input -3/-4/-5/-6 can search all possible three/four/five/six-center bonds" 310 write(*,*) "Input 0 can return to upper level menu" 311 read(*,"(a)") c1000tmp 312 313 if (c1000tmp(1:1)=='0') then 314 Return 315 else if (c1000tmp(1:1)/='-') then 316 call str2arr(c1000tmp,nbndcen,cenind) 317 318 do itime=1,ntime 319 if (wfntype==0.or.wfntype==3) then 320 PSmat=PSmattot 321 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 322 if (itime==1) PSmat=PSmattot 323 if (itime==2) PSmat=PSmata 324 if (itime==3) PSmat=PSmatb 325 end if 326 if (nbndcen>=8) write(*,*) "Please wait..." 327 call calcmultibndord(nbndcen,cenind,PSmat,nbasis,accum) !accum is pristine result without any factor 328 if (itime==1) bndordmix=accum 329 if (itime==2) bndordalpha=accum*2**(nbndcen-1) 330 if (itime==3) bndordbeta=accum*2**(nbndcen-1) 331 end do 332 if (wfntype==0.or.wfntype==3) then 333 write(*,"(a,f16.10)") " The multicenter bond order:",accum 334 !Normalized multicenter bond order, see Electronic Aromaticity Index for Large Rings DOI: 10.1039/C6CP00636A 335 !When it is negative, first obtain **(1/n) using its absolute value, then multiply it by -1 336 write(*,"(a,f16.10)") " The normalized multicenter bond order:",accum/abs(accum) * (abs(accum)**(1D0/nbndcen)) 337 338 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 339 write(*,"(a,f13.7)") " The multicenter bond order from alpha density matrix:",bndordalpha 340 write(*,"(a,f13.7)") " The multicenter bond order from beta density matrix: ",bndordbeta 341 totbndorder=bndordalpha+bndordbeta 342 write(*,"(a,f13.7)") " The sum of multicenter bond order from alpha and beta parts: ",totbndorder 343 write(*,"(a,f13.7)") " Above result in normalized form:",totbndorder/abs(totbndorder) * (abs(totbndorder)**(1D0/nbndcen)) 344 write(*,"(a,f13.7)") " The multicenter bond order from mixed alpha&beta density matrix:",bndordmix 345 write(*,"(a,f13.7)") " Above result in normalized form:",bndordmix/abs(bndordmix) * (abs(bndordmix)**(1D0/nbndcen)) 346 end if 347 348 else if (c1000tmp(1:1)=='-') then !Automatical search 349 read(c1000tmp,*) nbndcen 350 nbndcen=abs(nbndcen) 351 PSmat=PSmattot 352 !Search all combinations. Owing to simplicity and efficiency consideration, for open-shell system, compulsory to use mixed alpha&beta density matrix 353 if (wfntype==1.or.wfntype==2.or.wfntype==4) write(*,"(a)") "Note: The bond order considered here comes from mixed alpha&beta density matrix" 354 write(*,*) 355 write(*,*) "The bond order larger than which value should be outputted?" 356 write(*,*) "Input a value, e.g. 0.03" 357 read(*,*) thres 358 359 nfound=0 360 maxbndord=0D0 361 if (nbndcen/=3) write(*,*) "Note: The search may be not exhaustive. Please wait..." 362 if (nbndcen==3) then 363 !All combinations 364 do iatm=1,ncenter 365 do jatm=iatm+1,ncenter 366 do katm=jatm+1,ncenter 367 cenind(1)=iatm 368 cenind(2)=jatm 369 cenind(3)=katm 370 !clockwise and anticlockwise 371 do iseq=1,2 372 if (iseq==2) call invarr(cenind,1,nbndcen) 373 call calcmultibndord(nbndcen,cenind,PSmat,nbasis,accum) 374 if (accum>=thres) then 375 write(*,"(3i8,' Result:'f16.10)") cenind(1:nbndcen),accum 376 nfound=nfound+1 377 end if 378 if (accum>maxbndord) then 379 maxbndord=accum 380 maxcenind=cenind 381 end if 382 end do 383 end do 384 end do 385 end do 386 else if (nbndcen==4) then 387nthreads=getNThreads() 388!$OMP PARALLEL DO private(iatm,jatm,katm,latm,cenind,iseq,accum) shared(nfound) schedule(dynamic) NUM_THREADS(nthreads) 389 do iatm=1,ncenter 390 do jatm=iatm+1,ncenter 391 do katm=jatm+1,ncenter 392 do latm=katm+1,ncenter 393 cenind(1)=iatm 394 cenind(2)=jatm 395 cenind(3)=katm 396 cenind(4)=latm 397 do iseq=1,2 398 if (iseq==2) call invarr(cenind,1,nbndcen) 399 call calcmultibndord(nbndcen,cenind,PSmat,nbasis,accum) 400 if (accum>=thres) then 401 write(*,"(4i8,' Result:'f16.10)") cenind(1:nbndcen),accum 402 nfound=nfound+1 403 end if 404 if (accum>maxbndord) then 405 maxbndord=accum 406 maxcenind=cenind 407 end if 408 end do 409 end do 410 end do 411 end do 412 end do 413!$OMP end parallel do 414 else if (nbndcen==5) then 415 do iatm=1,ncenter 416nthreads=getNThreads() 417!$OMP PARALLEL DO private(jatm,katm,latm,matm,cenind,iseq,accum) shared(nfound) schedule(dynamic) NUM_THREADS(nthreads) 418 do jatm=iatm+1,ncenter 419 do katm=jatm+1,ncenter 420 do latm=katm+1,ncenter 421 do matm=latm+1,ncenter 422 cenind(1)=iatm 423 cenind(2)=jatm 424 cenind(3)=katm 425 cenind(4)=latm 426 cenind(5)=matm 427 do iseq=1,2 428 if (iseq==2) call invarr(cenind,1,nbndcen) 429 call calcmultibndord(nbndcen,cenind,PSmat,nbasis,accum) 430 if (accum>=thres) then 431 write(*,"(5i8,' Result:'f16.10)") cenind(1:nbndcen),accum 432 nfound=nfound+1 433 end if 434 if (accum>maxbndord) then 435 maxbndord=accum 436 maxcenind=cenind 437 end if 438 end do 439 end do 440 end do 441 end do 442 end do 443!$OMP end parallel do 444 end do 445 else if (nbndcen==6) then 446 do iatm=1,ncenter 447nthreads=getNThreads() 448!$OMP PARALLEL DO private(jatm,katm,latm,matm,cenind,iseq,accum) shared(nfound) schedule(dynamic) NUM_THREADS(nthreads) 449 do jatm=iatm+1,ncenter 450 do katm=jatm+1,ncenter 451 do latm=katm+1,ncenter 452 do matm=latm+1,ncenter 453 do natm=matm+1,ncenter 454 cenind(1)=iatm 455 cenind(2)=jatm 456 cenind(3)=katm 457 cenind(4)=latm 458 cenind(5)=matm 459 cenind(6)=natm 460 do iseq=1,2 461 if (iseq==2) call invarr(cenind,1,nbndcen) 462 call calcmultibndord(nbndcen,cenind,PSmat,nbasis,accum) 463 if (accum>=thres) then 464 write(*,"(6i8,' Result:'f16.10)") cenind(1:nbndcen),accum 465 nfound=nfound+1 466 end if 467 if (accum>maxbndord) then 468 maxbndord=accum 469 maxcenind=cenind 470 end if 471 end do 472 end do 473 end do 474 end do 475 end do 476 end do 477!$OMP end parallel do 478 write(*,"(a,i5,a,i5)") " Finished:",iatm,"/",ncenter 479 end do 480 end if 481 if (nfound==0) write(*,*) "No multi-center bonds above criteria were found" 482 write(*,*) 483 write(*,*) "The maximum bond order is" 484 if (nbndcen==3) write(*,"(3i8,' Result:'f16.10)") maxcenind(1:nbndcen),maxbndord 485 if (nbndcen==4) write(*,"(4i8,' Result:'f16.10)") maxcenind(1:nbndcen),maxbndord 486 if (nbndcen==5) write(*,"(5i8,' Result:'f16.10)") maxcenind(1:nbndcen),maxbndord 487 if (nbndcen==6) write(*,"(6i8,' Result:'f16.10)") maxcenind(1:nbndcen),maxbndord 488 end if 489end do 490end subroutine 491 492 493 494!!------ Bond order analysis in NAO basis 495!itask=1: Multicenter bond order calculation 496!itask=2: Decompose Wiberg bond order as atomic orbital pair and atomic shell pair contributions 497!NBO output file with DMNAO keyword should be used as input file 498subroutine bndordNAO(itask) 499use defvar 500use util 501implicit real*8 (a-h,o-z) 502integer cenind(12) 503integer,allocatable :: NAOinit(:),NAOend(:) 504character :: c80tmp*80,c80tmp2*80,c1000tmp*1000 505character,allocatable :: centername(:)*2 506real*8,allocatable :: DMNAO(:,:),DMNAOb(:,:),DMNAOtot(:,:) 507integer,allocatable :: NAOcen(:) 508character,allocatable :: NAOlang(:)*7,NAOcenname(:)*2,NAOtype(:)*3,NAOshtype(:)*2 509character*2 :: icenshname(100),jcenshname(100) !Record all shell type names in centers i and j, used for decomposition Wiberg bond order 510real*8,allocatable :: shcontri(:,:) 511 512open(10,file=filename,status="old") 513 514call loclabel(10,"NAO density matrix:",ifound,1) 515if (ifound==0) then 516 write(*,"(a)") "Error: Cannot found density matrix in NAO basis in the input file! Please read manual carefully" 517 write(*,*) 518else !Acquire number of NAOs and centers 519 call loclabel(10,"NATURAL POPULATIONS",ifound,1) 520 read(10,*) 521 read(10,*) 522 read(10,*) 523 read(10,*) 524 ilastspc=0 525 do while(.true.) !Find how many centers and how many NAOs. We need to carefully check where is ending 526 read(10,"(a)") c80tmp 527 if (c80tmp==' '.or.index(c80tmp,"low occupancy")/=0.or.index(c80tmp,"Population inversion found")/=0.or.index(c80tmp,"effective core potential")/=0) then 528 if (ilastspc==1) then 529 ncenter=iatm 530 numNAO=inao 531 exit 532 end if 533 ilastspc=1 !last line is space 534 else 535 read(c80tmp,*) inao,c80tmp2,iatm 536 ilastspc=0 537 end if 538 end do 539 write(*,"(' The number of atoms:',i10)") ncenter 540 write(*,"(' The number of NAOs: ',i10)") numNAO 541 allocate(NAOinit(ncenter),NAOend(ncenter),centername(ncenter),NAOcen(numNAO)) 542 !Get relationship between centers and NAO indices 543 call loclabel(10,"NATURAL POPULATIONS",ifound,1) 544 read(10,*) 545 read(10,*) 546 read(10,*) 547 read(10,*) 548 ilastspc=1 549 do while(.true.) 550 read(10,"(a)") c80tmp 551 if (c80tmp/=' ') then 552 read(c80tmp,*) inao,c80tmp2,iatm 553 NAOcen(inao)=iatm 554 if (ilastspc==1) NAOinit(iatm)=inao 555 ilastspc=0 556 else 557 NAOend(iatm)=inao 558 centername(iatm)=c80tmp2 559 if (iatm==ncenter) exit 560 ilastspc=1 561 end if 562 end do 563end if 564if (allocated(basstart)) deallocate(basstart,basend) 565allocate(basstart(ncenter),basend(ncenter)) 566basstart=NAOinit 567basend=NAOend 568 569if (itask==2) then !Load detailed NAO information 570 allocate(NAOcenname(numNAO),NAOtype(numNAO),NAOlang(numNAO),NAOshtype(numNAO)) 571 !Get relationship between center and NAO indices, as well as center name, then store these informationto NAOcen 572 !We delay to read NAO information, because in alpha and beta cases may be different 573 call loclabel(10,"NATURAL POPULATIONS",ifound,1) 574 read(10,*) 575 read(10,*) 576 read(10,*) 577 read(10,*) 578 do iatm=1,ncenter 579 do inao=NAOinit(iatm),NAOend(iatm) 580 read(10,"(a)") c80tmp 581 if (index(c80tmp,"Cor")/=0) then 582 NAOtype(inao)="Cor" 583 else if (index(c80tmp,"Val")/=0) then 584 NAOtype(inao)="Val" 585 else if (index(c80tmp,"Ryd")/=0) then 586 NAOtype(inao)="Ryd" 587 end if 588 read(c80tmp,*) c80tmp2,NAOcenname(inao),c80tmp2,NAOlang(inao),c80tmp2,c80tmp2 589 NAOshtype(inao)=c80tmp2(1:2) 590 end do 591 read(10,*) 592 end do 593end if 594 595!Read in density matrix in NAO basis. NAO information always use those generated by total density 596!DMNAO and AONAO is different for total, alpha and beta density 597allocate(DMNAO(numNAO,numNAO)) 598call loclabel(10,"NAO density matrix:",ifound) 599if (ifound==0) then 600 write(*,*) "Error: Cannot found NAO density matrix field in the input file" 601 write(*,*) 602 return 603end if 604!Check format before reading, NBO6 use different format to NBO3 605!For open-shell case, DMNAO doesn't print for total, so the first time loaded DMNAO is alpha(open-shell) or total(close-shell) 606read(10,"(a)") c80tmp 607backspace(10) 608nskipcol=16 !NBO3 609if (c80tmp(2:2)==" ") nskipcol=17 !NBO6 610call readmatgau(10,DMNAO,0,"f8.4 ",nskipcol,8,3) 611call loclabel(10," Alpha spin orbitals ",iopenshell) 612if (iopenshell==1) then 613 write(*,*) "This is an open-shell calculation" 614 allocate(DMNAOb(numNAO,numNAO),DMNAOtot(numNAO,numNAO)) 615 call loclabel(10,"******* Beta spin orbitals *******",ifound) 616 call loclabel(10,"NAO density matrix:",ifound,0) 617 call readmatgau(10,DMNAOb,0,"f8.4 ",nskipcol,8,3) 618 DMNAOtot=DMNAO+DMNAOb 619end if 620 621close(10) 622 623! call showmatgau(DMNAO,"Density matrix in NAO basis ",0,"f14.8",6) 624 625if (itask==1) then !Multi-center bond order calculation 626 do while(.true.) 627 write(*,*) 628 write(*,*) "Input atom indices, e.g. 3,4,7,8,10 (2~12 centers)" 629 write(*,*) "Input 0 can exit" 630 read(*,"(a)") c1000tmp 631 if (c1000tmp(1:1)=='0') then 632 deallocate(basstart,basend) 633 return 634 else 635 call str2arr(c1000tmp,nbndcen,cenind) 636 if (nbndcen>=7) write(*,*) "Please wait..." 637 638 if (iopenshell==0) then 639 call calcmultibndord(nbndcen,cenind,DMNAO,numNAO,bndord) 640 if (nbndcen==2) then 641 write(*,"(a,f16.10)") " The Wiberg bond order:",bndord 642 else 643 write(*,"(a,f16.10)") " The multicenter bond order:",bndord 644 write(*,"(a,f16.10)") " The normalized multicenter bond order:",bndord/abs(bndord) * (abs(bndord)**(1D0/nbndcen)) 645 end if 646 else if (iopenshell==1) then 647 call calcmultibndord(nbndcen,cenind,DMNAO,numNAO,bndordalpha) 648 call calcmultibndord(nbndcen,cenind,DMNAOb,numNAO,bndordbeta) 649 call calcmultibndord(nbndcen,cenind,DMNAOtot,numNAO,bndordmix) 650 bndordalpha=bndordalpha*2**(nbndcen-1) 651 bndordbeta=bndordbeta*2**(nbndcen-1) 652 totbndorder=bndordalpha+bndordbeta 653 if (nbndcen==2) then 654 write(*,"(a,f16.10)") " The bond order from alpha density matrix:",bndordalpha 655 write(*,"(a,f16.10)") " The bond order from beta density matrix: ",bndordbeta 656 write(*,"(a,f16.10)") " The sum of above two terms:",bndordalpha+bndordbeta 657 write(*,"(a,f16.10)") " The bond order from mixed alpha&beta density matrix: ",bndordmix 658 else 659 write(*,"(a,f13.7)") " The multicenter bond order from alpha density matrix:",bndordalpha 660 write(*,"(a,f13.7)") " The multicenter bond order from beta density matrix: ",bndordbeta 661 write(*,"(a,f13.7)") " The sum of multicenter bond order from alpha and beta parts: ",totbndorder 662 write(*,"(a,f13.7)") " Above result in normalized form:",totbndorder/abs(totbndorder) * (abs(totbndorder)**(1D0/nbndcen)) 663 write(*,"(a,f13.7)") " The multicenter bond order from mixed alpha&beta density matrix:",bndordmix 664 write(*,"(a,f13.7)") " Above result in normalized form:",bndordmix/abs(bndordmix) * (abs(bndordmix)**(1D0/nbndcen)) 665 end if 666 end if 667 end if 668 end do 669 670else if (itask==2) then !Decompose Wiberg bond order as atomic orbital pair contribution 671 if (iopenshell==1) then 672 write(*,*) "Error: This function currently only supports closed-shell system" 673 write(*,*) "Press ENTER to return" 674 read(*,*) 675 return 676 end if 677 write(*,"(a)") " Note: The threshold for printing contribution is controlled by ""bndordthres"" in settings.ini" 678 do while(.true.) 679 write(*,*) 680 write(*,*) "Input two atom indices, e.g. 3,4" 681 write(*,*) "Input 0 can exit" 682 read(*,"(a)") c1000tmp 683 if (c1000tmp(1:1)=='0') return 684 read(c1000tmp,*) iatm,jatm 685 !Construct name list of all shells for iatm and jatm 686 numicensh=1 687 icenshname(1)=NAOshtype(basstart(iatm)) 688 do ibas=basstart(iatm)+1,basend(iatm) 689 if (all(icenshname(1:numicensh)/=NAOshtype(ibas))) then 690 numicensh=numicensh+1 691 icenshname(numicensh)=NAOshtype(ibas) 692 end if 693 end do 694 numjcensh=1 695 jcenshname(1)=NAOshtype(basstart(jatm)) 696 do jbas=basstart(jatm)+1,basend(jatm) 697 if (all(jcenshname(1:numjcensh)/=NAOshtype(jbas))) then 698 numjcensh=numjcensh+1 699 jcenshname(numjcensh)=NAOshtype(jbas) 700 end if 701 end do 702! write(*,*) numicensh,numjcensh 703 allocate(shcontri(numicensh,numjcensh)) 704 shcontri=0 705 !Calculate Wiberg bond order and output worthnoting component 706 bndord=0 707 write(*,*) "Contribution from NAO pairs that larger than printing threshold:" 708 write(*,*) " Contri. NAO Center NAO type NAO Center NAO type" 709 do ibas=basstart(iatm),basend(iatm) 710 do ish=1,numicensh !Find the belonging shell index within this atom for NAO ibas 711 if (NAOshtype(ibas)==icenshname(ish)) exit 712 end do 713 do jbas=basstart(jatm),basend(jatm) 714 do jsh=1,numjcensh 715 if (NAOshtype(jbas)==jcenshname(jsh)) exit 716 end do 717 contri=DMNAO(ibas,jbas)**2 718 bndord=bndord+contri 719 shcontri(ish,jsh)=shcontri(ish,jsh)+contri 720 if (contri>bndordthres) write(*,"(f8.4,1x,i5,i5,'(',a,') ',a,'(',a,') ',a,'--- ',i5,i5,'(',a,') ',a,'(',a,') ',a)") contri,& 721 ibas,NAOcen(ibas),NAOcenname(ibas),NAOtype(ibas),NAOshtype(ibas),NAOlang(ibas),& 722 jbas,NAOcen(jbas),NAOcenname(jbas),NAOtype(jbas),NAOshtype(jbas),NAOlang(jbas) 723 end do 724 end do 725 write(*,*) 726 write(*,*) "Contribution from NAO shell pairs that larger than printing threshold:" 727 write(*,*) " Contri. Shell Center Type Shell Center Type" 728 do ish=1,numicensh 729 do jsh=1,numjcensh 730 if (shcontri(ish,jsh)>bndordthres) write(*,"(f8.4,1x,i5,i5,'(',a,') ',a,' --- ',i5,i5,'(',a,') ',a)") shcontri(ish,jsh),& 731 ish,NAOcen(basstart(iatm)),NAOcenname(basstart(iatm)),icenshname(ish),& 732 jsh,NAOcen(basstart(jatm)),NAOcenname(basstart(jatm)),jcenshname(jsh) 733 end do 734 end do 735 write(*,"(/,a,f8.4)") " Total Wiberg bond order:",bndord 736 deallocate(shcontri) 737 end do 738end if 739end subroutine 740 741 742!!---- A routine directly calculate multi-center bond order without complex things, shared by multicenter and bndordNAO 743!Note that for open-shell cases, the result varible "result" may then be multiplied by a proper factor 744subroutine calcmultibndord(nbndcen,cenind,PSmat,matdim,result) 745use defvar 746implicit real*8(a-h,o-z) 747real*8 PSmat(matdim,matdim) 748integer nbndcen,cenind(12) 749result=0D0 750if (nbndcen==2) then 751 do ib=basstart(cenind(2)),basend(cenind(2)) 752 do ia=basstart(cenind(1)),basend(cenind(1)) 753 result=result+PSmat(ia,ib)*PSmat(ib,ia) 754 end do 755 end do 756else if (nbndcen==3) then 757 do ic=basstart(cenind(3)),basend(cenind(3)) 758 do ib=basstart(cenind(2)),basend(cenind(2)) 759 do ia=basstart(cenind(1)),basend(cenind(1)) 760 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,ia) 761 end do 762 end do 763 end do 764else if (nbndcen==4) then 765 do id=basstart(cenind(4)),basend(cenind(4)) 766 do ic=basstart(cenind(3)),basend(cenind(3)) 767 do ib=basstart(cenind(2)),basend(cenind(2)) 768 do ia=basstart(cenind(1)),basend(cenind(1)) 769 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ia) 770 end do 771 end do 772 end do 773 end do 774else if (nbndcen==5) then 775 do ie=basstart(cenind(5)),basend(cenind(5)) 776 do id=basstart(cenind(4)),basend(cenind(4)) 777 do ic=basstart(cenind(3)),basend(cenind(3)) 778 do ib=basstart(cenind(2)),basend(cenind(2)) 779 do ia=basstart(cenind(1)),basend(cenind(1)) 780 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ie)*PSmat(ie,ia) 781 end do 782 end do 783 end do 784 end do 785 end do 786else if (nbndcen==6) then 787 do i_f=basstart(cenind(6)),basend(cenind(6)) 788 do ie=basstart(cenind(5)),basend(cenind(5)) 789 do id=basstart(cenind(4)),basend(cenind(4)) 790 do ic=basstart(cenind(3)),basend(cenind(3)) 791 do ib=basstart(cenind(2)),basend(cenind(2)) 792 do ia=basstart(cenind(1)),basend(cenind(1)) 793 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ie)*PSmat(ie,i_f)*PSmat(i_f,ia) 794 end do 795 end do 796 end do 797 end do 798 end do 799 end do 800else if (nbndcen==7) then 801 do i_g=basstart(cenind(7)),basend(cenind(7)) 802 do i_f=basstart(cenind(6)),basend(cenind(6)) 803 do ie=basstart(cenind(5)),basend(cenind(5)) 804 do id=basstart(cenind(4)),basend(cenind(4)) 805 do ic=basstart(cenind(3)),basend(cenind(3)) 806 do ib=basstart(cenind(2)),basend(cenind(2)) 807 do ia=basstart(cenind(1)),basend(cenind(1)) 808 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ie)*PSmat(ie,i_f)*PSmat(i_f,i_g)*PSmat(i_g,ia) 809 end do 810 end do 811 end do 812 end do 813 end do 814 end do 815 end do 816else if (nbndcen==8) then 817 do i_h=basstart(cenind(8)),basend(cenind(8)) 818 do i_g=basstart(cenind(7)),basend(cenind(7)) 819 do i_f=basstart(cenind(6)),basend(cenind(6)) 820 do ie=basstart(cenind(5)),basend(cenind(5)) 821 do id=basstart(cenind(4)),basend(cenind(4)) 822 do ic=basstart(cenind(3)),basend(cenind(3)) 823 do ib=basstart(cenind(2)),basend(cenind(2)) 824 do ia=basstart(cenind(1)),basend(cenind(1)) 825 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ie)*PSmat(ie,i_f)*PSmat(i_f,i_g)*PSmat(i_g,i_h)*PSmat(i_h,ia) 826 end do 827 end do 828 end do 829 end do 830 end do 831 end do 832 end do 833 end do 834else if (nbndcen==9) then 835 do i_i=basstart(cenind(9)),basend(cenind(9)) 836 do i_h=basstart(cenind(8)),basend(cenind(8)) 837 do i_g=basstart(cenind(7)),basend(cenind(7)) 838 do i_f=basstart(cenind(6)),basend(cenind(6)) 839 do ie=basstart(cenind(5)),basend(cenind(5)) 840 do id=basstart(cenind(4)),basend(cenind(4)) 841 do ic=basstart(cenind(3)),basend(cenind(3)) 842 do ib=basstart(cenind(2)),basend(cenind(2)) 843 do ia=basstart(cenind(1)),basend(cenind(1)) 844 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ie)*PSmat(ie,i_f)*PSmat(i_f,i_g)*PSmat(i_g,i_h)*PSmat(i_h,i_i)*PSmat(i_i,ia) 845 end do 846 end do 847 end do 848 end do 849 end do 850 end do 851 end do 852 end do 853 end do 854else if (nbndcen==10) then 855 itmp=0 856 ntot=basend(cenind(10))-basstart(cenind(10))+1 857 do i_j=basstart(cenind(10)),basend(cenind(10)) 858 itmp=itmp+1 859 write(*,"(' Finished:',i5,' /',i5)") itmp,ntot 860 do i_i=basstart(cenind(9)),basend(cenind(9)) 861 do i_h=basstart(cenind(8)),basend(cenind(8)) 862 do i_g=basstart(cenind(7)),basend(cenind(7)) 863 do i_f=basstart(cenind(6)),basend(cenind(6)) 864 do ie=basstart(cenind(5)),basend(cenind(5)) 865 do id=basstart(cenind(4)),basend(cenind(4)) 866 do ic=basstart(cenind(3)),basend(cenind(3)) 867 do ib=basstart(cenind(2)),basend(cenind(2)) 868 do ia=basstart(cenind(1)),basend(cenind(1)) 869 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ie)*PSmat(ie,i_f)*PSmat(i_f,i_g)*PSmat(i_g,i_h)*PSmat(i_h,i_i)*PSmat(i_i,i_j)*PSmat(i_j,ia) 870 end do 871 end do 872 end do 873 end do 874 end do 875 end do 876 end do 877 end do 878 end do 879 end do 880else if (nbndcen==11) then 881 itmp=0 882 ntot=( basend(cenind(11))-basstart(cenind(11))+1 ) * ( basend(cenind(10))-basstart(cenind(10))+1 ) 883 do i_k=basstart(cenind(11)),basend(cenind(11)) 884 do i_j=basstart(cenind(10)),basend(cenind(10)) 885 itmp=itmp+1 886 write(*,"(' Finished:',i5,' /',i5)") itmp,ntot 887 do i_i=basstart(cenind(9)),basend(cenind(9)) 888 do i_h=basstart(cenind(8)),basend(cenind(8)) 889 do i_g=basstart(cenind(7)),basend(cenind(7)) 890 do i_f=basstart(cenind(6)),basend(cenind(6)) 891 do ie=basstart(cenind(5)),basend(cenind(5)) 892 do id=basstart(cenind(4)),basend(cenind(4)) 893 do ic=basstart(cenind(3)),basend(cenind(3)) 894 do ib=basstart(cenind(2)),basend(cenind(2)) 895 do ia=basstart(cenind(1)),basend(cenind(1)) 896 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ie)*PSmat(ie,i_f)*PSmat(i_f,i_g)*PSmat(i_g,i_h)*PSmat(i_h,i_i)*PSmat(i_i,i_j)*PSmat(i_j,i_k)*PSmat(i_k,ia) 897 end do 898 end do 899 end do 900 end do 901 end do 902 end do 903 end do 904 end do 905 end do 906 end do 907 end do 908else if (nbndcen==12) then 909 itmp=0 910 ntot=( basend(cenind(12))-basstart(cenind(12))+1 ) * ( basend(cenind(11))-basstart(cenind(11))+1 ) * ( basend(cenind(10))-basstart(cenind(10))+1 ) 911 do i_l=basstart(cenind(12)),basend(cenind(12)) 912 do i_k=basstart(cenind(11)),basend(cenind(11)) 913 do i_j=basstart(cenind(10)),basend(cenind(10)) 914 itmp=itmp+1 915 write(*,"(' Finished:',i5,' /',i5)") itmp,ntot 916 do i_i=basstart(cenind(9)),basend(cenind(9)) 917 do i_h=basstart(cenind(8)),basend(cenind(8)) 918 do i_g=basstart(cenind(7)),basend(cenind(7)) 919 do i_f=basstart(cenind(6)),basend(cenind(6)) 920 do ie=basstart(cenind(5)),basend(cenind(5)) 921 do id=basstart(cenind(4)),basend(cenind(4)) 922 do ic=basstart(cenind(3)),basend(cenind(3)) 923 do ib=basstart(cenind(2)),basend(cenind(2)) 924 do ia=basstart(cenind(1)),basend(cenind(1)) 925 result=result+PSmat(ia,ib)*PSmat(ib,ic)*PSmat(ic,id)*PSmat(id,ie)*PSmat(ie,i_f)*PSmat(i_f,i_g)*PSmat(i_g,i_h)*PSmat(i_h,i_i)*PSmat(i_i,i_j)*PSmat(i_j,i_k)*PSmat(i_k,i_l)*PSmat(i_l,ia) 926 end do 927 end do 928 end do 929 end do 930 end do 931 end do 932 end do 933 end do 934 end do 935 end do 936 end do 937 end do 938end if 939end subroutine 940 941 942 943 944 945!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 946!!--------- Calculate Mulliken bond order 947subroutine mullikenbndord 948use defvar 949use util 950implicit real*8 (a-h,o-z) 951real*8 :: PSmata(nbasis,nbasis),PSmatb(nbasis,nbasis) 952real*8 :: bndmattot(ncenter,ncenter),bndmatb(ncenter,ncenter),bndmata(ncenter,ncenter) 953character selectyn 954bndmattot=0D0 955if (wfntype==0.or.wfntype==3) then 956 PSmata=Sbas*Ptot !Condensed to basis function matrix 957 do i=1,ncenter !Contract PSmata to Condensed to "Condensed to atoms" 958 do j=i+1,ncenter 959 bndmattot(i,j)=sum(PSmata(basstart(i):basend(i),basstart(j):basend(j))) 960 end do 961 end do 962 bndmattot=2*(bndmattot+transpose(bndmattot)) 963 forall (i=1:ncenter) bndmattot(i,i)=sum(bndmattot(i,:)) 964else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 965 bndmata=0D0 966 bndmatb=0D0 967 PSmata=Palpha*Sbas 968 PSmatb=Pbeta*Sbas 969 do i=1,ncenter 970 do j=i+1,ncenter 971 bndmata(i,j)=sum(PSmata(basstart(i):basend(i),basstart(j):basend(j))) 972 bndmatb(i,j)=sum(PSmatb(basstart(i):basend(i),basstart(j):basend(j))) 973 end do 974 end do 975 bndmata=2*(bndmata+transpose(bndmata)) 976 bndmatb=2*(bndmatb+transpose(bndmatb)) 977 forall (i=1:ncenter) bndmata(i,i)=sum(bndmata(i,:)) 978 forall (i=1:ncenter) bndmatb(i,i)=sum(bndmatb(i,:)) 979 bndmattot=bndmata+bndmatb 980end if 981 982write(*,"(' The absolute value of bond order >=',f10.6)") bndordthres 983itmp=0 984do i=1,ncenter 985 do j=i+1,ncenter 986 if (wfntype==0.or.wfntype==3) then 987 if (abs(bndmattot(i,j))>=bndordthres) then 988 itmp=itmp+1 989 write(*,"(' #',i5,':',5x,i5,a,i5,a,f14.8)") itmp,i,'('//a(i)%name//')',j,'('//a(j)%name//')',bndmattot(i,j) 990 end if 991 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 992 if (abs(bndmata(i,j)+bndmatb(i,j))>=bndordthres) then 993 itmp=itmp+1 994 write(*,"(' #',i5,':',i5,a,i5,a,' Alpha: ',f10.6,' Beta:',f10.6,' Total:',f10.6)") & 995 itmp,i,'('//a(i)%name//')',j,'('//a(j)%name//')',bndmata(i,j),bndmatb(i,j),bndmattot(i,j) 996 end if 997 end if 998 end do 999end do 1000write(*,*) 1001 1002!Between fragment 1003if (allocated(frag1)) then 1004 bndordfraga=0 1005 bndordfragb=0 1006 bndordfragtot=0 1007 do i=1,size(frag1) 1008 do j=1,size(frag2) 1009 bndordfraga=bndordfraga+bndmata(frag1(i),frag2(j)) 1010 bndordfragb=bndordfragb+bndmatb(frag1(i),frag2(j)) 1011 bndordfragtot=bndordfragtot+bndmattot(frag1(i),frag2(j)) 1012 end do 1013 end do 1014 if (wfntype==1.or.wfntype==2.or.wfntype==4) then 1015 write(*,"(' The Mulliken bond order between fragment 1 and 2:')") 1016 write(*,"(' Alpha:',f12.6,' Beta:',f12.6,' Total:',f12.6)") bndordfraga,bndordfragb,bndordfragtot 1017 else if (wfntype==0.or.wfntype==3) then 1018 write(*,"(' The Mulliken bond order between fragment 1 and 2:',f12.6)") bndordfragtot 1019 end if 1020 write(*,*) 1021end if 1022 1023write(*,*) "If output bond order matrix to bndmat.txt in current folder? (y/n)" 1024read(*,*) selectyn 1025if (selectyn=='y'.or.selectyn=='Y') then 1026 open(10,file="bndmat.txt",status="replace") 1027 write(10,*) "Note:The diagonal elements are the sum of corresponding row elements" 1028 if (wfntype==0.or.wfntype==3) then 1029 call showmatgau(bndmattot,"Mulliken bond order matrix",0,"f14.8",10) 1030 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 1031 call showmatgau(bndmata,"Mulliken bond order matrix for alpha electrons",0,"f14.8",10) 1032 call showmatgau(bndmatb,"Mulliken bond order matrix for beta electrons",0,"f14.8",10) 1033 call showmatgau(bndmattot,"Mulliken bond order matrix all electrons",0,"f14.8",10) 1034 end if 1035 close(10) 1036 write(*,*) "Result have been outputted to bndmat.txt in current folder" 1037 write(*,*) 1038end if 1039end subroutine 1040 1041 1042 1043!!--------- Decompose Mulliken bond order to MO contribution 1044subroutine decompMBO 1045use defvar 1046use util 1047implicit real*8 (a-h,o-z) 1048real*8,pointer :: ptmat(:,:) 1049 1050do while(.true.) 1051 write(*,*) "Input index of two atom (e.g. 3,5)" 1052 write(*,*) "Note: Input 0,0 can return to upper level menu" 1053 read(*,*) ind1,ind2 1054 1055 if (ind1==0.and.ind2==0) exit 1056 bndorda=0D0 1057 bndordb=0D0 1058 do itime=1,2 1059 if (itime==1) ptmat=>cobasa 1060 if (itime==2) ptmat=>cobasb 1061 if (itime==1.and.(wfntype==1.or.wfntype==4)) write(*,*) "Alpha orbitals:" 1062 if (itime==2.and.(wfntype==1.or.wfntype==4)) write(*,*) "Beta orbitals:" 1063 do imo=1,nbasis 1064 if (itime==1) irealmo=imo 1065 if (itime==2) irealmo=imo+nbasis 1066 if (MOocc(irealmo)==0D0) cycle 1067 accum=0D0 1068 do i=basstart(ind1),basend(ind1) 1069 do j=basstart(ind2),basend(ind2) 1070 accum=accum+MOocc(irealmo)*ptmat(i,imo)*ptmat(j,imo)*Sbas(i,j) 1071 end do 1072 end do 1073 if (itime==1) bndorda=bndorda+accum*2 1074 if (itime==2) bndordb=bndordb+accum*2 1075 write(*,"(' Orbital',i6,' Occ:',f10.6,' Energy:',f12.6,' contributes',f14.8)") imo,MOocc(irealmo),MOene(irealmo),accum*2 1076 end do 1077 if (wfntype==0.or.wfntype==2.or.wfntype==3) then 1078 write(*,"(' Total Mulliken bond order:',f14.8,/)") bndorda 1079 exit 1080 else if (wfntype==1.or.wfntype==4) then 1081 if (itime==1) write(*,"(' Mulliken bond order of all alpha electrons:',f14.8,/)") bndorda 1082 if (itime==2) write(*,"(' Mulliken bond order of all beta electrons:',f14.8,/)") bndordb 1083 if (itime==2) write(*,"(' Total Mulliken bond order:',f14.8,/)") bndorda+bndordb 1084 end if 1085 end do 1086end do 1087end subroutine 1088 1089 1090 1091!!----- Orbital occupancy-perturbed Mayer bond order (Decompose Mayer bond-order between two atoms to orbital contributions) 1092!--- J. Chem. Theory Comput. 2012, 8, 908, 914 1093!For simplicity, this routine only calculate Mayer bond for alpha and beta and then sum them up, don't concern mixed alpha+beta cases 1094subroutine OrbPertMayer 1095use defvar 1096implicit real*8 (a-h,o-z) 1097character orbtypechar*2 1098real*8 :: bndmata(ncenter,ncenter),bndmatb(ncenter,ncenter),bndmattot(ncenter,ncenter) 1099real*8,allocatable :: PSmattot(:,:),Ptottmp(:,:) 1100real*8,allocatable :: PSmata(:,:),PSmatb(:,:),Palphatmp(:,:),Pbetatmp(:,:) 1101bndmata=0D0 1102bndmatb=0D0 1103do while(.true.) 1104write(*,*) "Input indices of two atoms, e.g. 3,5" 1105 read(*,*) iatm,jatm 1106 if (iatm>=1.and.iatm<=ncenter.and.jatm>=1.and.jatm<=ncenter.and.iatm/=jatm) exit 1107 write(*,*) "Error: Invalid input, please input again" 1108end do 1109if (wfntype==0.or.wfntype==3) then !Close shell 1110 allocate(PSmattot(nbasis,nbasis),Ptottmp(nbasis,nbasis)) 1111 sumupvar=0D0 1112 do imo=0,nmo !Cycle all MOs 1113 Ptottmp=Ptot !Don't use Ptot to make troubles, because Ptot is a global array 1114 if (imo/=0) then !Calculate perturbed density. At the first time (imo=1), we don't pertube density matrix to yield original Mayer bond order 1115 if (MOocc(imo)<=1D-10) cycle 1116 do ibas=1,nbasis 1117 do jbas=1,nbasis 1118 Ptottmp(ibas,jbas)=Ptottmp(ibas,jbas)-MOocc(imo)*CObasa(ibas,imo)*CObasa(jbas,imo) 1119 end do 1120 end do 1121 end if 1122 PSmattot=matmul(Ptottmp,Sbas) !Calculate Mayer bond order based on Ptottmp 1123 bndordtot=0D0 1124 do ii=basstart(iatm),basend(iatm) 1125 do jj=basstart(jatm),basend(jatm) 1126 bndordtot=bndordtot+PSmattot(ii,jj)*PSmattot(jj,ii) 1127 end do 1128 end do 1129 if (imo==0) then 1130 beforepert=bndordtot 1131 write(*,"(' Mayer bond order before orbital occupancy-perturbation:',f12.6)") beforepert 1132 write(*,*) 1133 write(*,"(' Mayer bond order after orbital occupancy-perturbation:')") 1134 write(*,*) "Orbital Occ Energy Bond order Variance" 1135 else 1136 bndordvar=bndordtot-beforepert 1137 write(*,"(i6,f12.5,f11.5,2f12.6)") imo,MOocc(imo),MOene(imo),bndordtot,bndordvar 1138 sumupvar=sumupvar+bndordvar 1139 end if 1140 end do 1141 write(*,"(' Summing up occupancy perturbation from all orbitals:',f10.5)") sumupvar 1142 1143else if (wfntype==1.or.wfntype==2.or.wfntype==4) then !Open shell 1144 sumupvar=0D0 1145 allocate(PSmata(nbasis,nbasis),PSmatb(nbasis,nbasis),Palphatmp(nbasis,nbasis),Pbetatmp(nbasis,nbasis)) 1146 do imo=0,nmo 1147 Palphatmp=Palpha 1148 Pbetatmp=Pbeta 1149 if (imo/=0) then !The first time, we calculate actual Mayer bond order 1150 if (MOocc(imo)<=1D-10) cycle 1151 if (wfntype==1.or.wfntype==4) then 1152 if (imo<=nbasis) then !Alpha orbitals 1153 do ibas=1,nbasis 1154 do jbas=1,nbasis 1155 Palphatmp(ibas,jbas)=Palphatmp(ibas,jbas)-MOocc(imo)*CObasa(ibas,imo)*CObasa(jbas,imo) 1156 end do 1157 end do 1158 else !Beta orbitals, between nbasis+1 and nmo 1159 do ibas=1,nbasis 1160 do jbas=1,nbasis 1161 Pbetatmp(ibas,jbas)=Pbetatmp(ibas,jbas)-MOocc(imo)*CObasb(ibas,imo-nbasis)*CObasb(jbas,imo-nbasis) 1162 end do 1163 end do 1164 end if 1165 else if (wfntype==2) then !ROHF 1166 if (MOtype(imo)==0) then !Doubly occupied orbitals 1167 do ibas=1,nbasis 1168 do jbas=1,nbasis 1169 Palphatmp(ibas,jbas)=Palphatmp(ibas,jbas)-1D0*CObasa(ibas,imo)*CObasa(jbas,imo) 1170 Pbetatmp(ibas,jbas)=Pbetatmp(ibas,jbas)-1D0*CObasa(ibas,imo)*CObasa(jbas,imo) !For ROHF, Cobasb==Cobasa, and hence Cobasb is not allocated 1171 end do 1172 end do 1173 else if (MOtype(imo)==1) then !Alpha orbitals 1174 do ibas=1,nbasis 1175 do jbas=1,nbasis 1176 Palphatmp(ibas,jbas)=Palphatmp(ibas,jbas)-1D0*CObasa(ibas,imo)*CObasa(jbas,imo) 1177 end do 1178 end do 1179 end if 1180 end if 1181 end if 1182 PSmata=matmul(Palphatmp,Sbas) 1183 PSmatb=matmul(Pbetatmp,Sbas) 1184 bndorda=0D0 1185 bndordb=0D0 1186 do ii=basstart(iatm),basend(iatm) 1187 do jj=basstart(jatm),basend(jatm) 1188 bndorda=bndorda+PSmata(ii,jj)*PSmata(jj,ii) 1189 bndordb=bndordb+PSmatb(ii,jj)*PSmatb(jj,ii) 1190 end do 1191 end do 1192 bndorda=bndorda*2 1193 bndordb=bndordb*2 1194 if (imo==0) then 1195 beforepert=bndorda+bndordb 1196 write(*,"(' Mayer bond order before orbital occupancy-perturbation:')") 1197 write(*,"(' Alpha:',f12.6,' Beta:',f12.6,' Total:',f12.6)") bndorda,bndordb,bndorda+bndordb 1198 write(*,*) 1199 write(*,"(' Mayer bond order after orbital occupancy-perturbation:')") 1200 write(*,*) "Orbital Occ Energy Type Alpha Beta Total Variance" 1201 else 1202 bndordvar=bndorda+bndordb-beforepert 1203 if (MOtype(imo)==0) orbtypechar="AB" 1204 if (MOtype(imo)==1) orbtypechar="A " 1205 if (MOtype(imo)==2) orbtypechar="B " 1206 write(*,"(i6,f12.6,f11.5,2x,a,2x,3f10.5,3x,f10.5)") imo,MOocc(imo),MOene(imo),orbtypechar,bndorda,bndordb,bndorda+bndordb,bndordvar 1207 sumupvar=sumupvar+bndordvar 1208 end if 1209 end do 1210 write(*,"(' Summing up occupancy perturbation from all orbitals:',f10.5)") sumupvar 1211end if 1212write(*,*) 1213end subroutine 1214