1!!--------- Interface of orbital composition analysis 2subroutine compana 3use defvar 4implicit real*8 (a-h,o-z) 5!Automatically set proper radpot and sphpot for e.g. Becke, Hirshfeld, Hirshfeld-I 6radpotold=radpot 7sphpotold=sphpot 8if (iautointgrid==1) then 9 radpot=50 10 sphpot=170 11end if 12 13do while(.true.) 14 write(*,"(/,a,/)") " Citation of this module: Tian Lu, Feiwu Chen, Calculation of Molecular Orbital Composition, Acta Chim. Sinica, 69, 2393-2406" 15 write(*,*) " ================ Orbital composition analysis ===============" 16 write(*,*) "-10 Return" 17 if (allocated(CObasa)) then 18 write(*,*) "-2 Define fragment 2 (for option 4,5)" 19 write(*,*) "-1 Define fragment 1 (for option 1~6)" 20 write(*,*) "1 Orbital composition analysis with Mulliken partition" 21 write(*,*) "2 Orbital composition analysis with Stout-Politzer partition" 22 write(*,*) "3 Orbital composition analysis with Ros-Schuit (SCPA) partition" 23 write(*,*) "4 Print frag. 1 & inter-fragment compositions in all orbitals (Mulliken)" 24 write(*,*) "5 Print frag. 1 & inter-fragment compositions in all orbitals (Stout-Politzer)" 25 write(*,*) "6 Print frag. 1 compositions in all orbitals (SCPA)" 26 end if 27 write(*,*) "7 Orbital composition analysis by natural atomic orbital (NAO) method" 28 if (allocated(b)) write(*,*) "8 Calculate atom and fragment contributions by Hirshfeld method" 29 if (allocated(b)) write(*,*) "9 Calculate atom and fragment contributions by Becke method" 30 if (allocated(b)) write(*,*) "10 Calculate atom and fragment contributions by Hirshfeld-I method" 31 if (allocated(CObasa)) write(*,*) "100 Evaluate oxidation state by LOBA method" 32 read(*,*) icompana 33 34 if (icompana==-10) then 35 if (allocated(frag1)) deallocate(frag1) 36 if (allocated(frag2)) deallocate(frag2) 37 if (iautointgrid==1) then 38 radpot=radpotold 39 sphpot=sphpotold 40 end if 41 exit 42 else if (icompana==-1) then 43 call deffrag(1) 44 else if (icompana==-2) then 45 call deffrag(2) 46 else if (icompana==1.or.icompana==2.or.icompana==3) then 47 if (icompana==1) call orbcomponent(1) 48 if (icompana==2) call orbcomponent(2) 49 if (icompana==3) call orbcomponent(3) 50 else if (icompana==4.or.icompana==5.or.icompana==6) then 51 if (.not.allocated(frag1)) then 52 write(*,*) "Please use function -1 to define fragment #1" 53 write(*,*) 54 cycle 55 end if 56 if (.not.allocated(frag2).and.(icompana==4.or.icompana==5)) & 57 write(*,*) "Note: Fragment 2 was not be defined" 58 write(*,*) 59 if (icompana==4) call fragcomponent(1) 60 if (icompana==5) call fragcomponent(2) 61 if (icompana==6) call fragcomponent(3) 62 else if (icompana==7) then 63 call NAOcompana 64 else if (icompana==8) then 65 call spaceparorb(1) 66 else if (icompana==9) then 67 call spaceparorb(2) 68 else if (icompana==10) then 69 call Hirshfeld_I(2) 70 call spaceparorb(3) 71 else if (icompana==100) then 72 call LOBA 73 end if 74end do 75end subroutine 76 77 78 79!!-------------- Define fragment 1 & 2 and store to global array frag1 and frag2. "isel" denote define which fragment 80!If the fragment is empty, then will not allcoate frag1/2 when exit 81subroutine deffrag(isel) 82use defvar 83use util 84implicit real*8 (a-h,o-z) 85!fragtmp stores temporary basis index before exit setting interface, by default they are all zero, 86!if basis function 4,89,32 are in this fragment, then value 4,89,32 will be in uncertain three slots of fragtmp 87integer fragtmp(nbasis) 88integer termtmp(nbasis) !Store each time read serial number 89integer vectmp(nbasis) !used by "inv" command 90integer atmsellist(ncenter),bassellist(nbasis) 91character c200*200 92fragtmp=0 93if (isel==1.and.allocated(frag1)) then 94 fragtmp(1:size(frag1))=frag1(:) 95 deallocate(frag1) 96else if (isel==2.and.allocated(frag2)) then 97 fragtmp(1:size(frag2))=frag2(:) 98 deallocate(frag2) 99end if 10010 write(*,*) "Commands and examples:" 101write(*,*) "q: Save fragment and exit" 102write(*,*) "clean: Clean current fragment" 103write(*,*) "list: List basis functions in current fragment" 104write(*,*) "all: List all basis functions of current system" 105write(*,*) "addall: Add all basis functions to fragment" 106write(*,*) "inv: Invert selection, viz. delete existing basis functions add all other ones" 107write(*,*) "help: Print these help information again" 108write(*,*) "cond: Add basis functions satisfying given conditions to fragment" 109write(*,*) "a 2,5-8,12: Add all basis functions in atoms 2,5,6,7,8,12 to fragment" 110write(*,*) "s 2,5-8,12: Add all basis functions in shells 2,5,6,7,8,12 to fragment" 111write(*,*) "b 2,5-8,12: Add basis functions 2,5,6,7,8,12 to fragment" 112write(*,*) "da 5,7,11-13: Delete all basis functions in atoms 5,7,11,12,13 from fragment" 113write(*,*) "db 5,7,11-13: Delete basis functions 5,7,11,12,13 from fragment" 114 115do while(.true.) 116 read(*,"(a)") c200 117 if (c200(1:4)=="help") then 118 goto 10 119 else if (c200(1:6)=="addall") then 120 forall(i=1:nbasis) fragtmp(i)=i 121 write(*,*) "Done!" 122 else if (c200(1:1)=="q") then 123 numbas=count(fragtmp/=0) 124 if (isel==1.and.numbas>=1) allocate(frag1(numbas)) 125 if (isel==2.and.numbas>=1) allocate(frag2(numbas)) 126 write(*,*) "Basis function index in current fragment:" 127 if (numbas==0) then 128 write(*,*) "None" 129 else 130 do i=1,nbasis 131 if (fragtmp(i)/=0) write(*,"(i5)",advance="no") fragtmp(i) 132 end do 133 end if 134 write(*,*) 135 write(*,*) "Fragment is saved" 136 write(*,*) 137 if (numbas>=1) then 138 if (isel==1) frag1=pack(fragtmp,fragtmp/=0) 139 if (isel==2) frag2=pack(fragtmp,fragtmp/=0) 140 end if 141 exit 142 else if (c200(1:5)=="clean") then 143 fragtmp=0 144 write(*,*) "Done!" 145 else if (c200(1:3)=="all") then 146 write(*,*) "The basis functions with asterisk are those presented in current fragment" 147 do i=1,nbasis 148 if (any(fragtmp==i)) then 149 write(*,"('* Basis:',i6,' Shell:',i5,' Center:',i5,'(',a2,') Type: ',a)")& 150 i,basshell(i),bascen(i),a(bascen(i))%name,GTFtype2name(bastype(i)) 151 else 152 write(*,"(' Basis:',i6,' Shell:',i5,' Center:',i5,'(',a2,') Type: ',a)")& 153 i,basshell(i),bascen(i),a(bascen(i))%name,GTFtype2name(bastype(i)) 154 end if 155 end do 156 else if (c200(1:4)=="list") then 157 write(*,*) "Basis function in current fragment:" 158 if (all(fragtmp==0)) then 159 write(*,*) "None" 160 else 161 do i=1,nbasis 162 if (fragtmp(i)/=0) write(*,"(i5)",advance="no") fragtmp(i) 163 end do 164 write(*,*) 165 end if 166 else if (c200(1:3)=="inv") then 167 vectmp=fragtmp 168 fragtmp=0 169 itmp=0 170 do ibas=1,nbasis 171 if (all(vectmp/=ibas)) then 172 itmp=itmp+1 173 fragtmp(itmp)=ibas 174 end if 175 end do 176 write(*,*) "Done!" 177 else if (c200(1:4)=="cond") then 178 write(*,"(a)") " Note: You will be prompted to input three conditions in turn, & 179 the basis functions satisfying all conditions will be added to current fragment" 180 write(*,*) 181 write(*,*) "Condition 1: Input range of atoms, e.g. 2,5-8,12" 182 write(*,*) "To select all atoms, simply input ""a""" 183 read(*,"(a)") c200 184 if (index(c200,'a')/=0) then 185 nselatm=ncenter 186 forall(i=1:nselatm) atmsellist(i)=i 187 else 188 call str2arr(c200,nselatm,atmsellist) 189 if (any(atmsellist(1:nselatm)>ncenter)) then 190 write(*,*) "ERROR: One or more atom indices exceeded valid range!" 191 cycle 192 end if 193 end if 194 write(*,*) "Condition 2: Input range of basis functions, e.g. 2,5-8,12" 195 write(*,*) "To select all basis functions, simply input ""a""" 196 read(*,"(a)") c200 197 if (index(c200,'a')/=0) then 198 nselbas=nbasis 199 forall(i=1:nselbas) bassellist(i)=i 200 else 201 call str2arr(c200,nselbas,bassellist) 202 if (any(bassellist(1:nselbas)>nbasis)) then 203 write(*,*) "ERROR: One or more basis function indices exceeded valid range!" 204 cycle 205 end if 206 end if 207 write(*,*) "Condition 3: Choose basis function type, one of S,Y,Z,XY,YY,ZZZ..." 208 write(*,*) "You can also choose shell type, one of S, P, D, F, G, H" 209 write(*,*) """a"" means all types" 210 read(*,"(a)") c200 211 naddbas=0 212 do ibasidx=1,nselbas !Examine all basis functions 213 ibas=bassellist(ibasidx) 214 iadd=0 215 if ( any(fragtmp==ibas) ) cycle !Skip if already presented in current fragment 216 if ( all(atmsellist(1:nselatm)/=bascen(ibas)) ) cycle !Atom index condition 217 if ( index(c200,'a')==0) then !Basis function type condition 218 itype=bastype(ibas) 219 if ( trim(c200)=='P' .and. ( itype>=2.and.itype<=4 ) ) iadd=1 220 if ( trim(c200)=='D' .and. ( (itype>=-5 .and.itype<=-1 ).or.(itype>=5 .and.itype<=10) ) ) iadd=1 221 if ( trim(c200)=='F' .and. ( (itype>=-12.and.itype<=-6 ).or.(itype>=11.and.itype<=20) ) ) iadd=1 222 if ( trim(c200)=='G' .and. ( (itype>=-21.and.itype<=-13).or.(itype>=21.and.itype<=35) ) ) iadd=1 223 if ( trim(c200)=='H' .and. ( (itype>=-32.and.itype<=-22).or.(itype>=36.and.itype<=56) ) ) iadd=1 224 if ( trim(c200)==GTFtype2name(bastype(ibas)) ) iadd=1 !Inputted is detailed type 225 else 226 iadd=1 227 end if 228 if (iadd==1) then 229 do i=1,nbasis !Find space slot to save this basis function 230 if (fragtmp(i)==0) then 231 fragtmp(i)=ibas 232 naddbas=naddbas+1 233 exit 234 end if 235 end do 236 end if 237 end do 238 write(*,"(' Done!',i6,' new basis functions have been added to current fragment')") naddbas 239 240 else if (c200(1:2)=="a ".or.c200(1:2)=="s ".or.c200(1:2)=="b ".or.c200(1:2)=="db".or.c200(1:2)=="da") then 241 call str2arr(c200(3:),nterm,termtmp) 242 if (c200(1:2)=="a ") then 243 if (any(termtmp(1:nterm)<=0).or.any(termtmp(1:nterm)>ncenter)) then 244 write(*,*) "Atom index exceeded valid range! Ignoring..." 245 cycle 246 end if 247 do iatm=1,nterm 248 do ibas=basstart(termtmp(iatm)),basend(termtmp(iatm)) 249 if (all(fragtmp/=ibas)) then 250 do i=1,nbasis !Find an empty slot to record this basis function 251 if (fragtmp(i)==0) then 252 fragtmp(i)=ibas 253 exit 254 end if 255 end do 256 end if 257 end do 258 end do 259 else if (c200(1:2)=="s ") then 260 do ibas=1,nbasis 261 if (any(termtmp(1:nterm)==basshell(ibas)).and.all(fragtmp/=ibas)) then 262 do j=1,nbasis !Find an empty slot to record this basis function 263 if (fragtmp(j)==0) then 264 fragtmp(j)=ibas 265 exit 266 end if 267 end do 268 end if 269 end do 270 else if (c200(1:2)=="b ") then 271 do i=1,nterm 272 if (all(fragtmp/=termtmp(i)).and.termtmp(i)<=nbasis.and.termtmp(i)>0) then 273 do j=1,nbasis !Find empty slot to save this basis function 274 if (fragtmp(j)==0) then 275 fragtmp(j)=termtmp(i) 276 exit 277 end if 278 end do 279 end if 280 end do 281 else if (c200(1:2)=="db") then 282 do i=1,nterm 283 where(fragtmp==termtmp(i)) fragtmp=0 284 end do 285 else if (c200(1:2)=="da") then 286 do iatm=1,nterm 287 do ibas=basstart(termtmp(iatm)),basend(termtmp(iatm)) 288 where(fragtmp==ibas) fragtmp=0 289 end do 290 end do 291 end if 292 write(*,*) "Done!" 293 else 294 write(*,*) "Error: Unrecognized command, ignoring..." 295 end if 296end do 297if (allocated(frag1).and.allocated(frag2)) then 298 j=0 299 do i=1,size(frag2) 300 if (any(frag1==frag2(i))) then 301 write(*,"(' Warning: basis function',i6,' in fragment 2 also present in fragment 1!')") frag2(i) 302 j=1 303 end if 304 end do 305 if (j==1) then 306 write(*,*) "You must redefine fragment 1 or 2 to eliminate intersection set!" 307 write(*,*) 308 end if 309end if 310end subroutine 311 312 313 314!!--------- fragment and inter-fragment composition analysis 315!isel==1 fragment 1 & inter-fragment composition with Mulliken partition" 316!isel==2 fragment 1 & inter-fragment composition with Stout & Politzer partition" 317!isel==3 fragment 1 composition with Ros & Schuit (SCPA) partition" 318subroutine fragcomponent(isel) 319use defvar 320use util 321implicit real*8 (a-h,o-z) 322integer isel 323real*8,pointer :: tmpmat(:,:) 324real*8 :: ovpfrg12(nmo),ovpfrg12_1(nmo) 325character orbtype*2 326ovpfrg12=0D0 327ovpfrg12_1=0D0 328if (isel==1.or.isel==2) then 329 write(*,*) 330 write(*,"(a)") "Note:" 331 write(*,"(a)") """c^2"" means the square of coefficients of all basis functions within fragment 1" 332 write(*,"(a)") """Int.cross"" means cross term within fragment 1" 333 write(*,"(a,/)") """Ext.cross"" means the fragment 1 part of the total cross term between fragment 1 and all other basis functions" 334 write(*,"('Orb# Type Energy Occ c^2 Int.cross Ext.cross Total')") 335else if (isel==3) then 336 write(*,"('Orb# Type Energy Occ Composition')") 337end if 338 339do imo=1,nmo 340 if (imo<=nbasis) then 341 irealmo=imo 342 tmpmat=>cobasa 343 else if (imo>nbasis) then 344 irealmo=imo-nbasis 345 tmpmat=>cobasb 346 end if 347 floc=0D0 !Local term in fragment 1 348 crossint=0D0 349 crossext=0D0 !fragment 1 part of all overlap between fragment 1 and external atoms 350 if (isel==3) allsqr=sum(tmpmat(:,irealmo)**2) 351 do i=1,size(frag1) 352 ibas=frag1(i) 353 if (isel==1.or.isel==2) floc=floc+tmpmat(ibas,irealmo)**2 354 if (isel==3) floc=floc+tmpmat(ibas,irealmo)**2/allsqr 355 !Calculate cross term 356 if (isel==1.or.isel==2) then 357 do jbas=1,nbasis 358 if (jbas==ibas) cycle 359 crossij=tmpmat(ibas,irealmo)*tmpmat(jbas,irealmo)*Sbas(ibas,jbas) 360 if (any(frag1==jbas)) then 361 crossint=crossint+crossij !Internal cross term 362 else if (isel==1) then 363 crossext=crossext+crossij 364 else if (isel==2) then 365 tmpdenom=tmpmat(ibas,irealmo)**2+tmpmat(jbas,irealmo)**2 366 if (tmpdenom>1D-30) crossext=crossext+tmpmat(ibas,irealmo)**2/tmpdenom*crossij*2 367 end if 368 !Cross term between fragment 1 and 2. We assume there is no intersection set between frag1 and frag2 369 !Note: any(frag2==jbas) is a subset of .not.any(frag1==jbas), so ovpfrg12 is part of crossext 370 if (allocated(frag2).and.any(frag2==jbas)) then 371 ovpfrg12(imo)=ovpfrg12(imo)+2*crossij 372 if (isel==1) then 373 ovpfrg12_1(imo)=ovpfrg12_1(imo)+crossij 374 else if (isel==2) then 375 tmpdenom=tmpmat(ibas,irealmo)**2+tmpmat(jbas,irealmo)**2 376 if (tmpdenom>1D-30) ovpfrg12_1(imo)=ovpfrg12_1(imo)+tmpmat(ibas,irealmo)**2/tmpdenom*crossij*2 377 end if 378 end if 379 end do 380 end if 381 end do 382 if (MOtype(imo)==0) orbtype="AB" 383 if (MOtype(imo)==1) orbtype="A " 384 if (MOtype(imo)==2) orbtype="B " 385 if (isel==1.or.isel==2) write(*,"(i5,1x,a,f9.3,f7.3,f12.3,'%',f12.3,'%',f12.3,'%',f12.3,'%')") & 386 imo,orbtype,MOene(imo),MOocc(imo),floc*100,crossint*100,crossext*100,(floc+crossint+crossext)*100 387 if (isel==3) write(*,"(i5,1x,a,f9.3,f7.3,f18.6,'%')") imo,orbtype,MOene(imo),MOocc(imo),floc*100 388end do 389 390if (isel/=3.and.allocated(frag2)) then !Print cross term between fragment 1 and 2 391 write(*,*) 392 write(*,*) "Cross term between fragment 1 and 2 and their individual parts:" 393 write(*,"('Orb# Type Energy Occ Frag1 part Frag2 part Total')") 394 do imo=1,nmo 395 if (MOtype(imo)==0) orbtype="AB" 396 if (MOtype(imo)==1) orbtype="A " 397 if (MOtype(imo)==2) orbtype="B " 398 write(*,"(i5,1x,a,f9.3,f7.3,f14.4,'%',f14.4,'%',f14.4,'%')") & 399 imo,orbtype,MOene(imo),MOocc(imo),ovpfrg12_1(imo)*100,(ovpfrg12(imo)-ovpfrg12_1(imo))*100,ovpfrg12(imo)*100 400 end do 401end if 402end subroutine 403 404 405 406!!---------- Decompose orbitals to basis composition 407!isel=1 : Mulliken partition =2:Stout & Politzer partitio =3: Ros & Schuit (SCPA) partitionn" 408subroutine orbcomponent(isel) 409use defvar 410implicit real*8 (a-h,o-z) 411real*8 basloc(nbasis),bascross(nbasis),bastot(nbasis) 412real*8,pointer :: tmpmat(:,:) 413integer isel 414character orbtype*10 415do while(.true.) 416 write(*,*) "Input the orbital index to print composition, e.g. 4" 417 write(*,*) "Note: Input -1 can print basic information of all orbitals, input 0 to return" 418 read(*,*) ishowmo 419 if (ishowmo<-1.or.ishowmo>nmo) then 420 write(*,"('Orbital index should be in the range of 1 to',i6)") nmo 421 cycle 422 else if (ishowmo==0) then 423 exit 424 else if (ishowmo==-1) then !show all orbital information 425 do i=1,nmo 426 if (MOtype(i)==0) orbtype="Alpha&Beta" 427 if (MOtype(i)==1) orbtype="Alpha " 428 if (MOtype(i)==2) orbtype="Beta " 429 write(*,"('Orbital:',i5,' Energy(a.u.):',f14.8,' Occ:',f14.8,' Type: ',a)") i,MOene(i),MOocc(i),orbtype 430 end do 431 write(*,*) 432 else 433 write(*,*) 434 write(*,"('Threshold of absolute value: >',f12.6,'%')") compthres 435 if (MOtype(ishowmo)==0) orbtype="Alpha&Beta" 436 if (MOtype(ishowmo)==1) orbtype="Alpha " 437 if (MOtype(ishowmo)==2) orbtype="Beta " 438 if (ishowmo<=nbasis) tmpmat=>cobasa 439 if (ishowmo>nbasis) tmpmat=>cobasb 440 write(*,"('Orbital:',i5,' Energy(a.u.):',f14.8,' Occ:',f14.8,' Type: ',a)") ishowmo,MOene(ishowmo),MOocc(ishowmo),orbtype 441 if (isel==1.or.isel==2) write(*,"(' Basis Type Atom Shell Local Cross term Total ')") 442 if (isel==3) write(*,"(' Basis Type Atom Shell Composition')") 443 if (ishowmo>nbasis) ishowmo=ishowmo-nbasis !For wfntype==1.or.wfntype==4, change to #beta for cobasb 444 if (isel==3) allsqr=sum(tmpmat(:,ishowmo)**2) 445 bascross=0D0 446 do ibas=1,nbasis 447 basloc(ibas)=tmpmat(ibas,ishowmo)**2 448 if (isel==3) then 449 basloc(ibas)=basloc(ibas)/allsqr 450 else if (isel==1.or.isel==2) then 451 do jbas=1,nbasis 452 if (jbas==ibas) cycle 453 if (isel==1) then 454 bascross(ibas)=bascross(ibas)+tmpmat(ibas,ishowmo)*tmpmat(jbas,ishowmo)*Sbas(ibas,jbas) 455! write(*,*) ibas,jbas,tmpmat(ibas,ishowmo)*tmpmat(jbas,ishowmo)*Sbas(ibas,jbas) 456 else if (isel==2) then 457 tmp=tmpmat(ibas,ishowmo)**2+tmpmat(jbas,ishowmo)**2 458 if (tmp>1D-30) bascross(ibas)=bascross(ibas)+tmpmat(ibas,ishowmo)**2/tmp*2*tmpmat(ibas,ishowmo)*tmpmat(jbas,ishowmo)*Sbas(ibas,jbas) 459 end if 460 end do 461 end if 462 bastot(ibas)=basloc(ibas)+bascross(ibas) 463 if (abs(bastot(ibas))*100>compthres) then 464 if (isel==1.or.isel==2) write(*,"(i6,3x,a,i5,a,i5,f13.5,'%',f13.5,'%',f13.5,'%')") ibas,GTFtype2name(bastype(ibas)),bascen(ibas),'('//a(bascen(ibas))%name//')',& 465 basshell(ibas),basloc(ibas)*100,bascross(ibas)*100,bastot(ibas)*100 466 if (isel==3) write(*,"(i6,3x,a,i5,a,i5,f13.5,'%')") ibas,GTFtype2name(bastype(ibas)),bascen(ibas),'('//a(bascen(ibas))%name//')',& 467 basshell(ibas),bastot(ibas)*100 468 end if 469 end do 470 aboveloc=sum(basloc(:),abs(bastot)*100>compthres)*100 471 abovecross=sum(bascross(:),abs(bastot)*100>compthres)*100 472 if (isel==1.or.isel==2) write(*,"(' Sum up those listed above: ',f13.5,'%',f13.5,'%',f13.5,'%')") aboveloc,abovecross,aboveloc+abovecross 473 if (isel==1.or.isel==2) write(*,"(' Sum up all basis functions:',f13.5,'%',f13.5,'%',f13.5,'%')") sum(basloc(:))*100,sum(bascross(:))*100,sum(bastot(:))*100 474 if (isel==3) write(*,"(' Sum up those listed above: ',f13.5,'%')") sum(bastot(:),abs(bastot)*100>compthres)*100 475 if (isel==3) write(*,"(' Sum up all basis functions:',f13.5,'%')") sum(bastot(:))*100 476 write(*,*) 477 write(*,"(' Composition of each shell, threshold of absolute value: >',f12.6,'%')") compthres 478 do i=1,nshell 479 shellcom=0D0 480 do j=1,nbasis 481 if (basshell(j)==i) then 482 shellcom=shellcom+bastot(j) 483 iatm=bascen(j) 484 end if 485 end do 486 if (abs(shellcom)*100>compthres) write(*,"(' Shell',i6,' Type: ',a,' in atom',i5,'(',a,') :',f14.5,'%')") i,shtype2name(shtype(i)),iatm,a(iatm)%name,shellcom*100 487 end do 488 write(*,*) 489 write(*,*) "Composition of each atom:" 490 do i=1,ncenter 491! write(*,*) i,size(bastot),basstart(i),basend(i) 492 write(*,"(' Atom',i6,'(',a,') :',f12.6,'%')") i,a(i)%name,sum(bastot(basstart(i):basend(i)))*100 493 end do 494 if (allocated(frag1)) then 495 write(*,*) 496 write(*,"(' Composition of the fragment:',f12.6,'%')") sum(bastot(frag1(:)))*100 497 end if 498 write(*,*) 499 end if 500end do 501end subroutine 502 503 504 505 506!!!--- Partition orbital to atomic contribution by space partition methods 507!itype=1: Hirshfeld, itype=2: Becke, itype=3: Hirshfeld-I 508subroutine spaceparorb(itype) 509use defvar 510use util 511use function 512implicit real*8(a-h,o-z) 513type(content) gridorg(radpot*sphpot),gridatm(radpot*sphpot) 514real*8 resultvec(ncenter) 515real*8 allpotx(ncenter,radpot*sphpot),allpoty(ncenter,radpot*sphpot),allpotz(ncenter,radpot*sphpot),allpotw(ncenter,radpot*sphpot) 516real*8 tmpdens(radpot*sphpot),selfdens(radpot*sphpot),promol(radpot*sphpot),orbval(nmo),orbcomp(ncenter,nmo) 517integer,allocatable :: fragorbcomp(:) 518character orbtype*10,c2000tmp*2000 519real*8,external :: fdens_rad 520 521!Generate allpotx/y/z/w 522!allpotx(iatm,j) is x coordinate of the jth point for integrating center iatm 523write(*,*) "Initializing data, please wait... \(^_^)/" 524call gen1cintgrid(gridorg,iradcut) 525do iatm=1,ncenter 526 allpotx(iatm,:)=gridorg(:)%x+a(iatm)%x 527 allpoty(iatm,:)=gridorg(:)%y+a(iatm)%y 528 allpotz(iatm,:)=gridorg(:)%z+a(iatm)%z 529end do 530 531!allpotw combines Becke multi-center integration weight with Becke/Hirshfeld/Hirshfeld-I weight 532allpotw=0D0 533write(*,"(i7,' quadrature points are used for each atom')") radpot*sphpot 534call walltime(iwalltime1) 535if (itype==1.or.itype==3) then !Hirshfeld or Hirshfeld-I partition 536 write(*,*) 537 if (itype==1) then 538 write(*,*) "Hirshfeld analysis requests atomic densities, please select how to obtain them" 539 write(*,*) "1 Use build-in sphericalized atomic densities in free-states (recommended)" 540 write(*,"(a)") " 2 Provide wavefunction file of involved elements by yourself or invoke Gaussian to automatically calculate them" 541 read(*,*) ihirshmode 542 end if 543 if ((itype==1.and.ihirshmode==1).or.itype==3) then !Hirshfeld or Hirshfeld-I based on interpolation density 544 if (itype==1.and.ihirshmode==1) write(*,*) "Generating Hirshfeld atomic weighting functions at all grids..." 545 if (itype==3) write(*,*) "Generating Hirshfeld-I atomic weighting functions at all grids..." 546 do iatm=1,ncenter 547 promol=0D0 548 do jatm=1,ncenter 549 if (itype==1.and.ihirshmode==1) then !Hirshfeld based on interpolation of built-in atomic radius density 550nthreads=getNThreads() 551!$OMP parallel do shared(tmpdens) private(ipt) num_threads(nthreads) 552 do ipt=1+iradcut*sphpot,radpot*sphpot 553 tmpdens(ipt)=calcatmdens(jatm,allpotx(iatm,ipt),allpoty(iatm,ipt),allpotz(iatm,ipt),0) 554 end do 555!$OMP end parallel do 556 else if (itype==3) then !Hirshfeld-I based on refined atomic radial density 557nthreads=getNThreads() 558!$OMP parallel do shared(tmpdens) private(ipt) num_threads(nthreads) 559 do ipt=1+iradcut*sphpot,radpot*sphpot 560 tmpdens(ipt)=fdens_rad(jatm,allpotx(iatm,ipt),allpoty(iatm,ipt),allpotz(iatm,ipt)) 561 end do 562!$OMP end parallel do 563 end if 564 promol=promol+tmpdens 565 if (jatm==iatm) selfdens=tmpdens 566 end do 567 do i=1+iradcut*sphpot,radpot*sphpot !Get Hirshfeld weight of present atom 568 if (promol(i)/=0D0) then 569 allpotw(iatm,i)=selfdens(i)/promol(i) 570 else 571 allpotw(iatm,i)=0D0 572 end if 573 end do 574 allpotw(iatm,:)=allpotw(iatm,:)*gridorg(:)%value !Combine Hirshfeld weight with single-center integration weight 575 write(*,"(' Atom',i6,'(',a,') finished')") iatm,a(iatm)%name 576 end do 577 578 else if (itype==1.and.ihirshmode==2) then !Hirshfeld based on atomic .wfn file 579 call setpromol 580 do iatm=1,ncenter_org 581 promol=0D0 582 do jatm=1,ncenter_org 583 call dealloall 584 call readwfn(custommapname(jatm),1) 585nthreads=getNThreads() 586!$OMP parallel do shared(tmpdens) private(ipt) num_threads(nthreads) 587 do ipt=1+iradcut*sphpot,radpot*sphpot 588 tmpdens(ipt)=fdens(allpotx(iatm,ipt),allpoty(iatm,ipt),allpotz(iatm,ipt)) 589 end do 590!$OMP end parallel do 591 promol=promol+tmpdens 592 if (jatm==iatm) selfdens=tmpdens 593 end do 594 do i=1+iradcut*sphpot,radpot*sphpot !Get Hirshfeld weight of present atom 595 if (promol(i)/=0D0) then 596 allpotw(iatm,i)=selfdens(i)/promol(i) 597 else 598 allpotw(iatm,i)=0D0 599 end if 600 end do 601 allpotw(iatm,:)=allpotw(iatm,:)*gridorg(:)%value !Combine Hirshfeld weight with single-center integration weight 602 write(*,"(' Atom',i6,'(',a,') finished')") iatm,a_org(iatm)%name 603 end do 604 call dealloall 605 call readinfile(firstfilename,1) !Retrieve the firstly loaded file(whole molecule) in order to calculate real rho later 606 end if 607 608else if (itype==2) then !Becke partition 609 write(*,*) "Generating Becke weight..." 610 do iatm=1,ncenter !Cycle points of each atom 611 gridatm%x=gridorg%x+a(iatm)%x 612 gridatm%y=gridorg%y+a(iatm)%y 613 gridatm%z=gridorg%z+a(iatm)%z 614 call gen1cbeckewei(iatm,iradcut,gridatm,allpotw(iatm,:)) 615 allpotw(iatm,:)=allpotw(iatm,:)*gridorg%value !Combine Becke weight with single-center integration weight 616 write(*,"(' Atom',i6,'(',a,') finished')") iatm,a_org(iatm)%name 617 end do 618end if 619 620call walltime(iwalltime2) 621write(*,"(' Done! Initialization step took up wall clock time',i10,'s')") iwalltime2-iwalltime1 622write(*,*) 623 624do while(.true.) 625 write(*,*) "Now input the orbital index to print orbital composition, e.g. 5" 626 write(*,"(a)") " You can also input:" 627 if (.not.allocated(fragorbcomp)) then 628 write(*,"(a,i6)") " -9: Define fragment, current: undefined" 629 else 630 write(*,"(a,i6)") " -9: Redefine fragment, current number of atoms:",size(fragorbcomp) 631 end if 632 write(*,"(a)") " 0: Return" 633 write(*,"(a)") " -1: Print basic information of all orbitals" 634 write(*,"(a)") " -2: Print atom contribution to a range of orbitals" 635 write(*,"(a)") " -3: Print fragment contribution to a range of orbitals" 636 write(*,"(a)") " -4: Export composition of every atom in every orbital to orbcomp.txt" 637 read(*,*) ishowmo 638 if (ishowmo>nmo) then 639 write(*,"(' Error: Orbital index should within the range of 1 to',i6,/)") nmo 640 cycle 641 else if (ishowmo==0) then 642 exit 643 else if (ishowmo==-1) then !Show all orbital information 644 do i=1,nmo 645 if (MOtype(i)==0) orbtype="Alpha&Beta" 646 if (MOtype(i)==1) orbtype="Alpha " 647 if (MOtype(i)==2) orbtype="Beta " 648 write(*,"(' Orbital:',i5,' Energy(a.u.):',f14.8,' Occ:',f14.8,' Type: ',a)") i,MOene(i),MOocc(i),orbtype 649 end do 650 else if (ishowmo==-2.or.ishowmo==-3) then !Print atom/fragment contribution to specific range of orbitals 651 if ((.not.allocated(fragorbcomp)).and.ishowmo==-3) then 652 write(*,*) "Error: You must defined the fragment first!" 653 write(*,*) 654 cycle 655 end if 656 if (ishowmo==-2) then 657 write(*,*) "Input atom index, e.g. 4" 658 read(*,*) iatm 659 end if 660 write(*,*) "Input orbital range e.g. 20,25" 661 read(*,*) iorbbeg,iorbend 662 if (iorbbeg<=0) iorbbeg=1 663 if (iorbend>nmo) iorbend=nmo 664 write(*,"(' Orb# Type Ene(a.u.) Occ Composition Population')") 665 pop=0D0 666 do imo=iorbbeg,iorbend 667 if (MOtype(imo)==0) orbtype="Alpha&Beta" 668 if (MOtype(imo)==1) orbtype="Alpha " 669 if (MOtype(imo)==2) orbtype="Beta " 670 if (ishowmo==-2) then !Calculate only on atom 671 tmp=0D0 672nthreads=getNThreads() 673!$OMP parallel shared(tmp) private(ipot,value,tmpprivate) num_threads(nthreads) 674 tmpprivate=0D0 675!$OMP do schedule(dynamic) 676 do ipot=1+iradcut*sphpot,radpot*sphpot 677 if (allpotw(iatm,ipot)<1D-8) cycle !May lose 0.001% accuracy 678 value=fmo(allpotx(iatm,ipot),allpoty(iatm,ipot),allpotz(iatm,ipot),imo)**2 679 tmpprivate=tmpprivate+value*allpotw(iatm,ipot) 680 end do 681!$OMP end do 682!$OMP CRITICAL 683 tmp=tmp+tmpprivate 684!$OMP end CRITICAL 685!$OMP end parallel 686 else if (ishowmo==-3) then 687 tmp=0D0 688 do itmp=1,nfragorbcomp 689 iatm=fragorbcomp(itmp) 690nthreads=getNThreads() 691!$OMP parallel shared(tmp) private(ipot,value,tmpprivate) num_threads(nthreads) 692 tmpprivate=0D0 693!$OMP do schedule(dynamic) 694 do ipot=1+iradcut*sphpot,radpot*sphpot 695 if (allpotw(iatm,ipot)<1D-8) cycle !May lose 0.001% accuracy 696 value=fmo(allpotx(iatm,ipot),allpoty(iatm,ipot),allpotz(iatm,ipot),imo)**2 697 tmpprivate=tmpprivate+value*allpotw(iatm,ipot) 698 end do 699!$OMP end do 700!$OMP CRITICAL 701 tmp=tmp+tmpprivate 702!$OMP end CRITICAL 703!$OMP end parallel 704 end do 705 end if 706 write(*,"(i5,1x,a,f13.4,f9.3,f11.3,'%',f15.6)") imo,orbtype,MOene(imo),MOocc(imo),tmp*100,MOocc(imo)*tmp 707 pop=pop+MOocc(imo)*tmp 708 end do 709 if (iorbend-iorbbeg>0) write(*,"(a,f12.6)") " Population of this atom in these orbitals:",pop 710 711 else if (ishowmo==-4) then !Export composition of every atom in every orbital to orbcomp.txt in current folder 712 write(*,*) "Calculating, please wait..." 713 orbcomp=0 714nthreads=getNThreads() 715!$OMP parallel do shared(orbcomp) private(iatm,ipot,orbval) num_threads(nthreads) schedule(dynamic) 716 do iatm=1,ncenter 717 do ipot=1+iradcut*sphpot,radpot*sphpot 718 if (allpotw(iatm,ipot)<1D-9) cycle !May lose 0.001% accuracy 719 call orbderv(1,1,nmo,allpotx(iatm,ipot),allpoty(iatm,ipot),allpotz(iatm,ipot),orbval) 720 orbcomp(iatm,:)=orbcomp(iatm,:)+orbval(:)**2*allpotw(iatm,ipot) 721 end do 722 end do 723!$OMP end parallel do 724 open(10,file="orbcomp.txt",status="replace") 725 do imo=1,nmo 726 write(10,"(' Orbital',i6)") imo 727 do iatm=1,ncenter 728 write(10,"(i6,f12.6,' %')") iatm,orbcomp(iatm,imo)*100 729 end do 730 end do 731 close(10) 732 write(*,*) "Done! orbcomp.txt has been exported to current folder." 733 734 else if (ishowmo==-9) then !Define fragment 735 if (allocated(fragorbcomp)) then 736 write(*,*) "Atoms in current fragment:" 737 write(*,"(13i6)") fragorbcomp 738 write(*,"(a)") " Input 0 to keep unchanged, or redefine fragment, e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment" 739 else 740 write(*,"(a)") " Input atomic indices to define fragment. e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment" 741 end if 742 read(*,"(a)") c2000tmp 743 if (c2000tmp(1:1)=='0') then 744 continue 745 else if (c2000tmp(1:1)==" ") then 746 deallocate(fragorbcomp) 747 else 748 if (allocated(fragorbcomp)) deallocate(fragorbcomp) 749 call str2arr(c2000tmp,nfragorbcomp) 750 allocate(fragorbcomp(nfragorbcomp)) 751 call str2arr(c2000tmp,nfragorbcomp,fragorbcomp) 752 end if 753 754 else 755 if (MOtype(ishowmo)==0) orbtype="Alpha&Beta" 756 if (MOtype(ishowmo)==1) orbtype="Alpha " 757 if (MOtype(ishowmo)==2) orbtype="Beta " 758 write(*,"(' Orbital:',i5,' Energy(a.u.):',f14.8,' Occ:',f14.8,' Type: ',a)") ishowmo,MOene(ishowmo),MOocc(ishowmo),orbtype 759 write(*,*) "Please wait..." 760 accum=0D0 761 do iatm=1,ncenter 762 tmp=0D0 763nthreads=getNThreads() 764!$OMP parallel shared(tmp) private(ipot,value,tmpprivate) num_threads(nthreads) 765 tmpprivate=0 766!$OMP do schedule(dynamic) 767 do ipot=1+iradcut*sphpot,radpot*sphpot 768 if (allpotw(iatm,ipot)<1D-8) cycle !May lose 0.001% accuracy 769 value=fmo(allpotx(iatm,ipot),allpoty(iatm,ipot),allpotz(iatm,ipot),ishowmo)**2 770 tmpprivate=tmpprivate+value*allpotw(iatm,ipot) 771 end do 772!$OMP end do 773!$OMP CRITICAL 774 tmp=tmp+tmpprivate 775!$OMP end CRITICAL 776!$OMP end parallel 777 accum=accum+tmp 778 resultvec(iatm)=tmp 779 end do 780 write(*,"(' The sum of contributions before normalization',11x,f12.6,'%',/)") accum*100 781 if (allocated(fragorbcomp)) then 782 fragresult=sum(resultvec(fragorbcomp)) 783 end if 784 write(*,*) "Contributions after normalization:" 785 do iatm=1,ncenter 786 write(*,"(' Atom',i6,'(',a,') :',f12.6,'%')") iatm,a(iatm)%name,resultvec(iatm)/accum*100 787 end do 788 if (allocated(fragorbcomp)) write(*,"(/,' Fragment contribution:',f12.6,'%')") fragresult/accum*100 789 end if 790 write(*,*) 791end do 792end subroutine 793 794 795 796 797 798!!------------ Orbital composition analysis by NAO method 799subroutine NAOcompana 800use defvar 801use util 802implicit real*8 (a-h,o-z) 803integer :: ioutmode=1,ifragend=0,iorboutcomp !The orbital index selected 804integer :: ispinmode=0 !0/1=close/open-shell 805integer numNAO 806integer,allocatable :: fragidx(:),NAOinit(:),NAOend(:),NAOcen(:),termtmp(:) 807character :: c80tmp*80,c80tmp2*80,c200*200 !type2char(0:2)=(/"Cor","Val","Ryd"/) 808character,allocatable :: NAOlang(:)*7,NAOcenname(:)*2,NAOtype(:)*3,centername(:)*2,NAOshtype(:)*2 809real*8 :: outcrit=0.5D0 810real*8,allocatable :: NAOMO(:,:) 811 812open(10,file=filename,status="old") 813nmo=0 814call loclabel(10,"NBsUse=",ifound) !Gaussian may eliminate some linear dependency basis functions, so MO may be smaller than numNAO. NBsUse always equals to actual number of MO 815if (ifound==1) then 816 read(10,"(a)") c80tmp 817 !For DKH case, G09 may output such as RelInt: Using uncontracted basis, NUniq= 180 NBsUse= 180 0.100E-13, this nbsuse is meaningless, use next nbsuse 818 if (index(c80tmp,"RelInt")/=0) call loclabel(10,"NBsUse=",ifound) 819 itmp=index(c80tmp,'=') 820 read(c80tmp(itmp+1:),*) nmo 821end if 822call loclabel(10,"MOs in the NAO basis:",ifound,1) 823if (ifound==0) then 824 write(*,"(a)") " Error: Cannot found MOs in NAO basis in the input file, the input file is invalid for this function! Please read manual carefully" 825 write(*,*) 826 return 827else !Acquire number of NAOs and centers 828 call loclabel(10,"NATURAL POPULATIONS",ifound,1) 829 read(10,*) 830 read(10,*) 831 read(10,*) 832 read(10,*) 833 ilastspc=0 834 do while(.true.) !Find how many center and how many NAOs. We need to carefully check where is ending 835 read(10,"(a)") c80tmp 836 if (c80tmp==''.or.index(c80tmp,"low occupancy")/=0.or.index(c80tmp,"Population inversion found")/=0.or.index(c80tmp,"effective core potential")/=0) then 837 if (ilastspc==1) then 838 ncenter=iatm 839 numNAO=inao 840 exit 841 end if 842 ilastspc=1 !last line is space 843 else 844 read(c80tmp,*) inao,c80tmp2,iatm 845 ilastspc=0 846 end if 847 end do 848 if (nmo==0) nmo=numNAO 849 write(*,"(' The number of atoms:',i10)") ncenter 850 write(*,"(' The number of NAOs: ',i10)") numNAO 851 write(*,"(' The number of MOs : ',i10)") nmo 852 allocate(fragidx(numNAO),termtmp(numNAO)) 853 ifragend=0 !ifragend is acutal ending position of fragidx 854 allocate(NAOinit(ncenter),NAOend(ncenter),centername(ncenter)) 855 allocate(NAOcen(numNAO),NAOcenname(numNAO),NAOtype(numNAO),NAOlang(numNAO),NAOMO(numNAO,nmo),NAOshtype(numNAO)) 856 !Get relationship between center and NAO indices, as well as center name, then store these informationto NAOcen 857 !We delay to read NAO information, because in alpha and beta cases may be different 858 call loclabel(10,"NATURAL POPULATIONS",ifound,1) 859 read(10,*) 860 read(10,*) 861 read(10,*) 862 read(10,*) 863 ilastspc=1 864 do while(.true.) 865 read(10,"(a)") c80tmp 866 if (c80tmp/='') then 867 read(c80tmp,*) inao,c80tmp2,iatm 868 NAOcen(inao)=iatm 869 if (ilastspc==1) NAOinit(iatm)=inao 870 ilastspc=0 871 else 872 NAOend(iatm)=inao 873 centername(iatm)=c80tmp2 874 if (iatm==ncenter) exit 875 ilastspc=1 876 end if 877 end do 878end if 879 880!Determine if this file is close or open shell calculation 881call loclabel(10,"******* Alpha spin orbitals *******",ispinmode,1) 882if (ispinmode==0) write(*,*) "This is close-shell calculation" 883if (ispinmode==1) write(*,*) "This is open-shell calculation" 884 885!Interface 886NAOcompmaincyc: do while(.true.) 887do while(.true.) 888 write(*,*) 889 write(*,*) "-10 Return" 890 write(*,*) "-1 Define fragment (for option 1)" 891 write(*,*) "0 Show composition of an orbital" 892 write(*,*) "1 Show fragment contribution in a range of orbitals" 893 if (ioutmode==0) write(*,"(a)") " 2 Select output mode (for option 0), current: Show all NAOs" 894 if (ioutmode==1) write(*,"(a)") " 2 Select output mode (for option 0), current: Only show core and valence NAOs" 895 if (ioutmode==2) write(*,"(a,f6.2,'%')") " 2 Select output mode (for option 0), current: Show NAOs whose contributions are >",outcrit 896 if (ispinmode==1) write(*,*) "3 Switch spin type, current: Alpha" 897 if (ispinmode==2) write(*,*) "3 Switch spin type, current: Beta" 898 read(*,*) isel 899 900 if (isel==-10) then 901 close(10) 902 return 903 else if (isel==-1) then 90420 write(*,*) "Commands and examples:" 905 write(*,*) """q"" : Save setting and exit" 906 write(*,*) """help"" : Print these help content again" 907 write(*,*) """clean"": Clean current fragment" 908 write(*,*) """list"" : List NAOs in current fragment" 909 write(*,*) """all"" : Print all NAOs of current system" 910 write(*,*) """addall"": Add all NAOs to fragment" 911 write(*,*) """a 2,5-8,12"": Add all NAOs in atoms 2,5,6,7,8,12 to fragment" 912 write(*,*) """b 2,5-8,12"": Add NAOs 2,5,6,7,8,12 to fragment" 913 write(*,*) """da 5,7,11-13"": Delete all NAOs in atoms 5,7,11,12,13 from fragment" 914 write(*,*) """db 5,7,11-13"": Delete NAOs 5,7,11,12,13 from fragment" 915 do while(.true.) 916 read(*,"(a)") c200 917 if (c200(1:4)=="help") then 918 goto 20 919 else if (c200(1:6)=="addall") then 920 forall(i=1:numNAO) fragidx(i)=i 921 ifragend=numNAO 922 write(*,*) "Done!" 923 else if (c200(1:1)=="q") then 924 write(*,*) 925 write(*,*) "Fragment is saved, indices of NAOs in current fragment:" 926 if (ifragend==0) then 927 write(*,*) "None" 928 else 929 write(*,"(12i6)") fragidx(1:ifragend) 930 end if 931 exit 932 else if (c200(1:5)=="clean") then 933 fragidx=0 934 ifragend=0 935 write(*,*) "Done!" 936 else if (c200(1:3)=="all") then 937 if (ispinmode==0) then 938 call loclabel(10,"NATURAL POPULATIONS",ifound,1) 939 else if (ispinmode==1) then 940 call loclabel(10,"NATURAL POPULATIONS",ifound,1) !Total density 941 read(10,*) 942 call loclabel(10,"NATURAL POPULATIONS",ifound,0) !Alpha density 943 else if (ispinmode==2) then 944 call loclabel(10,"NATURAL POPULATIONS",ifound,1) !Total density 945 read(10,*) 946 call loclabel(10,"NATURAL POPULATIONS",ifound,0) !Alpha density 947 read(10,*) 948 call loclabel(10,"NATURAL POPULATIONS",ifound,0) !Beta density 949 end if 950 read(10,*) 951 read(10,*) 952 read(10,*) 953 read(10,*) 954 write(*,*) "Note: The ones with asterisks are those presented in current fragment" 955 write(*,*) " NAO# Atom Label Type Occupancy Energy" 956 itmp=0 957 do iatm=1,ncenter 958 do i=NAOinit(iatm),NAOend(iatm) 959 itmp=itmp+1 960 read(10,"(a)") c200 961 if (any(fragidx(1:ifragend)==itmp)) then 962 write(*,"(' *',a)") trim(c200) 963 else 964 write(*,"(' ',a)") trim(c200) 965 end if 966 end do 967 read(10,*) 968 end do 969 else if (c200(1:4)=="list") then 970 write(*,*) "NAO indices in current fragment:" 971 if (ifragend==0) then 972 write(*,*) "None" 973 else 974 write(*,"(12i6)") fragidx(1:ifragend) 975 end if 976 977 else if (c200(1:2)=="a ".or.c200(1:2)=="s ".or.c200(1:2)=="b ".or.c200(1:2)=="db".or.c200(1:2)=="da") then 978 call str2arr(c200(3:),nterm,termtmp) 979 if (c200(1:2)=="a ".or.c200(1:2)=="da") then !Check sanity of the input 980 if (any(termtmp(1:nterm)<=0).or.any(termtmp(1:nterm)>ncenter)) then 981 write(*,*) "ERROR: Atom index exceeded valid range! Ignoring..." 982 cycle 983 end if 984 else if (c200(1:2)=="b ".or.c200(1:2)=="db") then 985 if (any(termtmp(1:nterm)<=0).or.any(termtmp(1:nterm)>numNAO)) then 986 write(*,*) "ERROR: NAO index exceeded valid range! Ignoring..." 987 cycle 988 end if 989 end if 990 !Modify fragidx according to the selection 991 if (c200(1:2)=="a ") then 992 do iterm=1,nterm 993 iatm=termtmp(iterm) 994 do ibas=NAOinit(iatm),NAOend(iatm) 995 if ( any(fragidx(1:ifragend)==ibas) ) cycle 996 ifragend=ifragend+1 997 fragidx(ifragend)=ibas 998 end do 999 end do 1000 else if (c200(1:2)=="b ") then 1001 do ibas=1,nterm 1002 if ( any(fragidx(1:ifragend)==termtmp(ibas)) ) cycle 1003 ifragend=ifragend+1 1004 fragidx(ifragend)=termtmp(ibas) 1005 end do 1006 else if (c200(1:2)=="db") then 1007 do itmp=1,nterm 1008 do iscan=1,ifragend 1009 if (fragidx(iscan)==termtmp(itmp)) then 1010 fragidx(iscan:ifragend-1)=fragidx(iscan+1:ifragend) 1011 ifragend=ifragend-1 1012 exit 1013 end if 1014 if (iscan==ifragend) write(*,"('Note: NAO with index of',i7,' does not present in current fragment')") termtmp(itmp) 1015 end do 1016 end do 1017 else if (c200(1:2)=="da") then 1018 do itmp=1,nterm 1019 iatm=termtmp(itmp) 1020 do iNAOincen=NAOinit(iatm),NAOend(iatm) !Cycle NAO index in itmp atom 1021 do iscan=1,ifragend !Scan fragidx array to find out iNAOincen 1022 if (fragidx(iscan)==iNAOincen) then 1023 fragidx(iscan:ifragend-1)=fragidx(iscan+1:ifragend) 1024 ifragend=ifragend-1 1025 exit 1026 end if 1027 if (iscan==ifragend) write(*,"('Note: NAO with index of',i7,' does not present in current fragment')") iNAOincen 1028 end do 1029 end do 1030 end do 1031 end if 1032 write(*,*) "Done!" 1033 else 1034 write(*,*) "Error: Cannot recognize the command you inputted" 1035 end if 1036 1037 end do !End of fragment definition module 1038 1039 else if (isel==0) then 1040 exit 1041 1042 else if (isel==1) then 1043 if (ifragend==0) then 1044 write(*,*) "Error: You haven't defined fragment or the fragment is empty!" 1045 else 1046 write(*,*) "Input orbital range to be outputted e.g. 1,10" 1047 write(*,"(a,i7)") " Note: Should within 1 to",nmo 1048 read(*,*) iorblow,iorbhigh 1049 if (iorbhigh<iorblow.or.iorblow<=0.or.iorbhigh>nmo) then 1050 write(*,*) "Error: The range you inputted is invalid!" 1051 cycle 1052 else 1053 exit 1054 end if 1055 end if 1056 1057 else if (isel==2) then 1058 write(*,*) "0 Show all NAOs" 1059 write(*,*) "1 Only show core and valence NAOs" 1060 write(*,*) "2 Show NAOs whose contribution is larger than specified criteria" 1061 read(*,*) ioutmode 1062 if (ioutmode==2) then 1063 write(*,*) "Input criteria in percentage, e.g. 0.5" 1064 read(*,*) outcrit 1065 end if 1066 1067 else if (isel==3) then 1068 if (ispinmode==1) then 1069 ispinmode=2 1070 else 1071 ispinmode=1 1072 end if 1073 end if 1074end do 1075 1076!Start analysis! 1077if (isel==0.or.isel==1) then 1078 !Before analysis, read in NAO information according to spin mode 1079 if (ispinmode==0) then 1080 call loclabel(10,"NATURAL POPULATIONS",ifound,1) 1081 else if (ispinmode==1) then 1082 call loclabel(10,"NATURAL POPULATIONS",ifound,1) !Total density 1083 read(10,*) 1084 call loclabel(10,"NATURAL POPULATIONS",ifound,0) !Alpha density 1085 else if (ispinmode==2) then 1086 call loclabel(10,"NATURAL POPULATIONS",ifound,1) !Total density 1087 read(10,*) 1088 call loclabel(10,"NATURAL POPULATIONS",ifound,0) !Alpha density 1089 read(10,*) 1090 call loclabel(10,"NATURAL POPULATIONS",ifound,0) !Beta density 1091 end if 1092 read(10,*) 1093 read(10,*) 1094 read(10,*) 1095 read(10,*) 1096 do iatm=1,ncenter 1097 do inao=NAOinit(iatm),NAOend(iatm) 1098 read(10,"(a)") c80tmp 1099 if (index(c80tmp,"Cor")/=0) then 1100 NAOtype(inao)="Cor" 1101 else if (index(c80tmp,"Val")/=0) then 1102 NAOtype(inao)="Val" 1103 else if (index(c80tmp,"Ryd")/=0) then 1104 NAOtype(inao)="Ryd" 1105 end if 1106 read(c80tmp,*) c80tmp2,NAOcenname(inao),c80tmp2,NAOlang(inao),c80tmp2,c80tmp2 1107 NAOshtype(inao)=c80tmp2(1:2) 1108 end do 1109 read(10,*) 1110 end do 1111 !Read MO in NAO basis 1112 call loclabel(10,"MOs in the NAO basis:",ifound,0) !Don't rewind 1113 !Check format before reading, NBO6 use different format to NBO3 1114 read(10,"(a)") c80tmp 1115 backspace(10) 1116 if (c80tmp(2:2)==" ") then !NBO6 1117 call readmatgau(10,NAOMO,0,"f8.4 ",17,8,3) 1118 else !NBO3 1119 call readmatgau(10,NAOMO,0,"f8.4 ",16,8,3) 1120 end if 1121end if 1122 1123if (isel==0) then 1124 do while(.true.) 1125 write(*,*) "Analyze which orbital? (Input 0 can return)" 1126 read(*,*) iorboutcomp 1127 if (iorboutcomp==0) exit 1128 if (iorboutcomp<0.or.iorboutcomp>nmo) then 1129 write(*,"(a,i7)") "Error: The orbital index should between 1 and",nmo 1130 cycle 1131 end if 1132 1133 write(*,*) 1134 if (ioutmode==2) write(*,"(a,f6.2,a)") "Note: All NAOs whose contribution <=",outcrit,"% are ignored" 1135 if (ispinmode==1) write(*,*) "Below are composition of alpha orbitals" 1136 if (ispinmode==2) write(*,*) "Below are composition of beta orbitals" 1137 write(*,*) " NAO# Center Label Type Composition" 1138 sumcomp=0D0 1139 do inao=1,numNAO 1140 tmpcomp=NAOMO(inao,iorboutcomp)**2*100D0 1141 if (ioutmode==1.and.NAOtype(inao)=="Ryd") cycle !skip Ryd 1142 if (ioutmode==2.and.tmpcomp<=outcrit) cycle 1143 write(*,"( i8,i5,'(',a,')',4x,a,2x,a,'(',a,')',f14.6,'%' )") inao,NAOcen(inao),NAOcenname(inao),NAOlang(inao),NAOtype(inao),NAOshtype(inao),tmpcomp 1144 sumcomp=sumcomp+tmpcomp 1145 end do 1146 write(*,"(' Summing up the compositions listed above:',f14.6,'%')") sumcomp 1147 if (ioutmode==1) write(*,"( ' Rydberg composition:',f14.6,'%')") 100D0-sumcomp 1148 write(*,*) 1149 write(*,*) "Condensed above result to atoms:" 1150 write(*,*) " Center Composition" 1151 do icen=1,ncenter 1152 sumcomp=0D0 1153 do inao=NAOinit(icen),NAOend(icen) 1154 tmpcomp=NAOMO(inao,iorboutcomp)**2*100D0 1155 if (ioutmode==1.and.NAOtype(inao)=="Ryd") cycle !skip Ryd 1156 if (ioutmode==2.and.tmpcomp<=outcrit) cycle 1157 sumcomp=sumcomp+tmpcomp 1158 end do 1159 write(*,"( i6,'(',a,')',f12.6,'%' )") icen,centername(icen),sumcomp 1160 end do 1161 write(*,*) 1162 end do 1163else if (isel==1) then 1164 if (ispinmode==1) write(*,*) "Below are composition of alpha orbitals" 1165 if (ispinmode==2) write(*,*) "Below are composition of beta orbitals" 1166 write(*,*) " Orb.# Core Valence Rydberg Total" 1167 do imo=iorblow,iorbhigh 1168 sumcompcor=0D0 1169 sumcompval=0D0 1170 sumcompryd=0D0 1171 do itmp=1,ifragend 1172 inao=fragidx(itmp) 1173 tmpcomp=NAOMO(inao,imo)**2*100D0 1174 if (NAOtype(inao)=="Cor") sumcompcor=sumcompcor+tmpcomp 1175 if (NAOtype(inao)=="Val") sumcompval=sumcompval+tmpcomp 1176 if (NAOtype(inao)=="Ryd") sumcompryd=sumcompryd+tmpcomp 1177 end do 1178 sumcomptot=sumcompcor+sumcompval+sumcompryd 1179 write(*,"(i6,5x,4(f12.6,'%'))") imo,sumcompcor,sumcompval,sumcompryd,sumcomptot 1180 end do 1181end if 1182 1183end do NAOcompmaincyc 1184 1185end subroutine 1186 1187 1188 1189!!---------- Localized orbital bonding analysis (LOBA) based on SCPA partition 1190!Input file should contains localized MO 1191subroutine LOBA 1192use defvar 1193use util 1194implicit real*8 (a-h,o-z) 1195real*8 oxdstat(ncenter) 1196integer,allocatable :: fragLOBA(:) 1197character c2000tmp*2000 1198write(*,*) "Citation: Phys. Chem. Chem. Phys., 11, 11297-11304 (2009)" 1199write(*,*) "Note: SCPA method is employed here for calculating orbital composition" 1200write(*,*) 1201nfragLOBA=0 1202do while(.true.) 1203 oxdstat=a%charge 1204 write(*,*) "Input percentage threshold to determine oxidation state (50 is commonly used)" 1205 write(*,*) "Input -1 can define fragment" 1206 write(*,*) "Input 0 can exit" 1207 read(*,*) thres 1208 if (thres==0) then 1209 return 1210 else if (thres==-1) then 1211 write(*,*) "Input the atom indices constituting the fragment" 1212 write(*,*) "e.g. 1,3-5,9" 1213 read(*,"(a)") c2000tmp 1214 if (allocated(fragLOBA)) deallocate(fragLOBA) 1215 call str2arr(c2000tmp,nfragLOBA) 1216 allocate(fragLOBA(nfragLOBA)) 1217 call str2arr(c2000tmp,nfragLOBA,fragLOBA) 1218 write(*,*) "Done!" 1219 write(*,*) 1220 cycle 1221 else 1222 thres=thres/100 1223 do iatm=1,ncenter 1224 do imo=1,nmo 1225 if (MOtype(imo)==0.or.MOtype(imo)==1) then 1226 tmpnorm=sum(cobasa(:,imo)**2) 1227 compval=sum(cobasa(basstart(iatm):basend(iatm),imo)**2)/tmpnorm 1228 else if (MOtype(imo)==2) then 1229 tmpnorm=sum(cobasb(:,imo-nbasis)**2) 1230 compval=sum(cobasb(basstart(iatm):basend(iatm),imo-nbasis)**2)/tmpnorm 1231 end if 1232 if (compval>thres) oxdstat(iatm)=oxdstat(iatm)-MOocc(imo) 1233! if (MOocc(imo)>0.and.compval>thres) write(*,"(' Orbital',i4,' belongs to atom',i4)") imo,iatm 1234 end do 1235 write(*,"(' Oxidation state of atom',i4,'(',a,') :',i3)") iatm,a(iatm)%name,nint(oxdstat(iatm)) 1236 end do 1237 write(*,"(' The sum of oxidation states:',i4)") sum(nint(oxdstat)) 1238 if (allocated(fragLOBA)) then 1239 oxdfrag=0 1240 do iatmidx=1,nfragLOBA 1241 iatm=fragLOBA(iatmidx) 1242 oxdfrag=oxdfrag+a(iatm)%charge 1243 end do 1244 do imo=1,nmo 1245 compval=0 1246 if (MOtype(imo)==0.or.MOtype(imo)==1) then 1247 tmpnorm=sum(cobasa(:,imo)**2) 1248 else if (MOtype(imo)==2) then 1249 tmpnorm=sum(cobasb(:,imo-nbasis)**2) 1250 end if 1251 do iatmidx=1,nfragLOBA 1252 iatm=fragLOBA(iatmidx) 1253 if (MOtype(imo)==0.or.MOtype(imo)==1) then 1254 compval=compval+sum(cobasa(basstart(iatm):basend(iatm),imo)**2)/tmpnorm 1255 else if (MOtype(imo)==2) then 1256 compval=compval+sum(cobasb(basstart(iatm):basend(iatm),imo-nbasis)**2)/tmpnorm 1257 end if 1258 end do 1259 if (compval>thres) oxdfrag=oxdfrag-MOocc(imo) 1260 end do 1261 write(*,"(' Oxidation state of the fragment:',i4)") nint(oxdfrag) 1262 end if 1263 write(*,*) 1264 end if 1265end do 1266 1267end subroutine 1268