1!!------- Analyze or visualize hole-electron distribution, transition dipole moment and transition density 2!ROHF,RODFT are assumed to be impossible to be ground state 3!itype=1 is normal mode 4!itype=2 is specifically used to calculate delta_r 5!itype=3 is specifically used to generate NTOs 6!itype=4 is specifically used to Calculate inter-fragment charger transfer 7subroutine hetransdipdens(itype) 8use defvar 9use util 10use function 11implicit real*8 (a-h,o-z) 12integer :: itype,idomag=0 13logical,allocatable :: skippair(:) 14integer,allocatable :: excdir(:) !excdir=1 means ->, =2 means <- 15integer,allocatable :: orbleft(:),orbright(:) !denote the actual MO at the left/right side in the excitation data (1:nbasis=alpha/spatial, nbasis+1:2*nbasis=beta) 16real*8,allocatable :: exccoeff(:),exccoeffbackup(:) !Coefficient of an orbital pair transition. exccoeffbackup is used to backup, because users can modify the coefficients 17real*8,allocatable :: holegrid(:,:,:),elegrid(:,:,:),holeeleovlp(:,:,:),transdens(:,:,:),holecross(:,:,:),elecross(:,:,:),Cele(:,:,:),Chole(:,:,:),magtrdens(:,:,:,:) 18real*8,allocatable :: dipcontri(:,:) !(1/2/3,iexc) contribution of orbital pairs "iexc" to transition dipole moment in X/Y/Z 19character c200tmp*200,c80tmp*80,transmodestr*200,leftstr*80,rightstr*80,strtmp1*10,strtmp2*10,strtmp3*10,selectyn 20character,save :: excitfilename*200=" " 21integer walltime1,walltime2 22real*8 time_end,time_endtmp,orbval(nmo),wfnderv(3,nmo) 23real*8,allocatable :: tmparr1(:),tmparr2(:) !Arrays for temporary use 24real*8,allocatable :: tdmata(:,:),tdmatb(:,:) !Transition density matrix in basis function representation of alpha/total and beta electrons 25real*8,allocatable :: bastrdip(:,:) !Contribution from each basis function to transition dipole moment, the first index 1/2/3=X/Y/Z 26real*8,allocatable :: trdipmatbas(:,:,:),trdipmatatm(:,:,:) !Transition dipole moment matrix in basis function / atom representation, the first index 1/2/3=X/Y/Z 27real*8,allocatable :: bastrpopa(:),bastrpopb(:) !Transition population of each basis function 28!Store dipole moment integral between all GTFs, the matrix is symmetry hence stored as 1D-array to save memory. The first index is 1/2/3=X/Y/Z 29real*8,allocatable :: GTFdipint(:,:) 30!Store information of another excitation 31integer,allocatable :: orbleft2(:),orbright2(:),excdir2(:) 32real*8,allocatable :: exccoeff2(:) 33!Other 34real*8,allocatable :: cubx(:),cuby(:),cubz(:) 35! real*8 eigvecmat(nbasis,nbasis),eigvalarr(nbasis) 36 37if (.not.allocated(CObasa)) then 38 write(*,*) "Error: The input file does not contain basis function information!" 39 write(*,*) "Please read manual to make clear which kinds of input file could be used!" 40 write(*,*) "Press ENTER to exit" 41 read(*,*) 42 return 43end if 44 45if (excitfilename==" ") then 46 write(*,"(a)") " Input the path of the Gaussian/ORCA output file or plain text file containing excitation data, e.g. c:\a.out" 47 do while(.true.) 48 read(*,"(a)") excitfilename 49 inquire(file=excitfilename,exist=alive) 50 if (alive) exit 51 write(*,*) "Cannot find this file, input again" 52 end do 53else 54 write(*,"(' Loading ',a)") trim(excitfilename) 55end if 56open(10,file=excitfilename,status="old") 57call loclabel(10,"Gaussian",igauout,1,50) 58call loclabel(10,"O R C A",iORCAout,1,50) 59 60if (igauout==1) then !Gaussian output file 61 write(*,*) "Analyzing the file..." 62 call loclabel(10,"Excitation energies and oscillator strengths:") 63 read(10,*) 64 nstates=0 !The number of transition modes 65 do while(.true.) 66 call loclabel(10,"Excited State",ifound,0) 67 if (ifound==1) then 68 nstates=nstates+1 69 read(10,*) 70 else 71 exit 72 end if 73 end do 74 write(*,"(' There are',i5,' transition modes, analyze which one? e.g. 2')") nstates 75 do while(.true.) 76 read(*,*) iexcitmode 77 if (iexcitmode>0.and.iexcitmode<=nstates) exit 78 write(*,*) "Error: The index exceeded valid range, input again" 79 end do 80 call loclabel(10,"Excitation energies and oscillator strengths:") 81 do itmp=1,iexcitmode !Move to the iexcitmode field 82 call loclabel(10,"Excited State",ifound,0) 83 read(10,*) 84 end do 85 backspace(10) 86 read(10,"(a)") transmodestr 87 iexcmulti=1 !Multiplicity of the excited state 88 if (index(transmodestr,"Triplet")/=0) iexcmulti=3 89 do i=10,70 90 if (transmodestr(i:i+1)=="eV") exit 91 end do 92 read(transmodestr(i-10:i-1),*) excenergy 93else if (iORCAout==1) then !ORCA output file 94 write(*,*) "Analyzing the file..." 95 call loclabel(10,"Number of roots to be determined",ifound) 96 if (ifound==0) then 97 write(*,*) "Error: It seems that this is not a output file of CIS/TDA/TD task" 98 write(*,*) "Press Enter to exit" 99 read(*,*) 100 return 101 end if 102 read(10,"(50x,i7)") nstates 103 write(*,"(' There are',i5,' transition modes, analyze which one? e.g. 2')") nstates 104 do while(.true.) 105 read(*,*) iexcitmode 106 if (iexcitmode>0.and.iexcitmode<=nstates) exit 107 write(*,*) "Error: The index exceeded valid range, input again" 108 end do 109 iexcmulti=1 !Multiplicity of the excited state 110 if (wfntype==0.or.wfntype==3) then 111 call loclabel(10,"Generation of triplets") 112 read(10,"(a)") c80tmp 113 if (index(c80tmp," on ")/=0) then 114 write(*,*) "Load which kind of excited states?" 115 write(*,*) "1: Singlet 3: Triplet" 116 read(*,*) iexcmulti 117 end if 118 end if 119 write(*,*) "Loading data, please wait..." 120 call loclabel(10,"the weight of the individual excitations are printed") 121 if (iexcmulti==3) then !When triplets=on, ORCA calculate both singlet and triplet excited state, now move to the latter 122 read(10,*) 123 call loclabel(10,"the weight of the individual excitations are printed",ifound,0) 124 end if 125 do itmp=1,iexcitmode !Move to the iexcitmode field 126 call loclabel(10,"STATE",ifound,0) 127 read(10,*) 128 end do 129 backspace(10) 130 read(10,"(a)") transmodestr 131 do i=10,70 132 if (transmodestr(i:i+1)=="eV") exit 133 end do 134 read(transmodestr(i-10:i-1),*) excenergy 135else 136 rewind(10) 137 read(10,*) iexcmulti,excenergy !The first line should be multiplicity of excited state and excitation energy 138end if 139write(*,"(' Multiplicity of this excited state is',i4)") iexcmulti 140write(*,"(' Excitation energy is',f12.7,' eV')") excenergy 141 142!Count how many orbital pairs are involved in this transition mode 143nexcitorb=0 144do while(.true.) 145 read(10,"(a)",iostat=ierror) c80tmp 146 if ((index(c80tmp,'<-')==0.and.index(c80tmp,'->')==0).or.ierror/=0) exit 147 nexcitorb=nexcitorb+1 148end do 149allocate(excdir(nexcitorb),orbleft(nexcitorb),orbright(nexcitorb),exccoeff(nexcitorb),exccoeffbackup(nexcitorb)) 150if (igauout==1.or.iORCAout==1) then !Rewind to head of the entry 151 call loclabel(10,trim(transmodestr)) 152 read(10,*) 153else 154 rewind(10) 155 read(10,*) 156end if 157write(*,"(a,i8,a)") " There are",nexcitorb," orbital pairs in this transition mode" 158 159!Load transition detail. Notice that for unrestricted case, A and B are separately recorded in input file, & 160!however after loading, they are combined as single index, namely if orbital index is larger than nbasis, then it is B, else A 161if (iORCAout==1) then 162 !Worthnotingly, in at least ORCA 4.0, de-excitation is not separately outputted as <-, but combined into -> 163 !Here we still determine <-, because hopefully Neese may change the convention of ORCA output in the future... 164 do itmp=1,nexcitorb 165 read(10,"(a)") c80tmp 166 excdir(itmp)=-1 167 if (index(c80tmp,'->')/=0) excdir(itmp)=1 !means -> 168 if (index(c80tmp,'<-')/=0) excdir(itmp)=2 !means <- 169 if (excdir(itmp)==-1) then 170 write(*,"(a)") " Error: This output file does not correspond to a CIS or TD or TDA task, configuration information cannot be found!" 171 write(*,*) "Press Enter to exit" 172 read(*,*) 173 return 174 end if 175 do isign=1,80 !Find position of <- or -> 176 if (c80tmp(isign:isign)=='-'.or.c80tmp(isign:isign)=='<') exit 177 end do 178 !Process left side of <- or -> 179 read(c80tmp(:isign-1),"(a)") leftstr 180 read(leftstr(:len_trim(leftstr)-1),*) orbleft(itmp) 181 orbleft(itmp)=orbleft(itmp)+1 !ORCA counts orbital from 0 rather than 1!!! 182 if (index(leftstr,'b')/=0) orbleft(itmp)=orbleft(itmp)+nbasis 183 !Process right side of <- or -> 184 read(c80tmp(isign+2:),*) rightstr 185 read(rightstr(:len_trim(rightstr)-1),*) orbright(itmp) 186 orbright(itmp)=orbright(itmp)+1 187 if (index(rightstr,'b')/=0) orbright(itmp)=orbright(itmp)+nbasis 188 iTDA=index(c80tmp,'c=') 189 if (iTDA/=0) then !CIS, TDA task, configuration coefficients are presented 190 read(c80tmp(iTDA+2:iTDA+13),*) exccoeff(itmp) 191 else !TD task, configuration coefficients are not presented. Contribution of i->a and i<-a are summed up and outputted as i->a 192 if (itmp==1) then 193 write(*,"(a)") " Warning: For TD task, ORCA does not print configuration coefficients but only print corresponding contributions of each orbital pair, & 194 in this case Multiwfn determines configuration coefficients simply as square root of contribution values. However, this treatment is & 195 evidently inappropriate and the result is nonsense when de-excitation is significant (In this situation you have to use TDA-DFT instead)" 196 write(*,*) "If you really want to proceed, press ENTER to continue" 197 read(*,*) 198 end if 199 read(c80tmp(23:32),*) tmpval 200 if (tmpval<0) excdir(itmp)=2 !Negative contribution is assumed to be de-excitation (of course this is not strict since -> and <- have been combined together) 201 exccoeff(itmp)=dsqrt(abs(tmpval)) 202 end if 203 !Although for closed-shell ground state, ORCA still outputs coefficients as normalization to 100%, & 204 !However, in order to follow the Gaussian convention, we change the coefficient as normalization to 50% 205 if (wfntype==0.or.wfntype==3) exccoeff(itmp)=exccoeff(itmp)/dsqrt(2D0) 206 end do 207else !Gaussian output file or plain text file 208 do itmp=1,nexcitorb 209 read(10,"(a)") c80tmp 210 excdir(itmp)=1 !means -> 211 if (index(c80tmp,'<-')/=0) excdir(itmp)=2 !means <- 212 do isign=1,80 !Find position of <- or -> 213 if (c80tmp(isign:isign)=='-'.or.c80tmp(isign:isign)=='<') exit 214 end do 215 !Process left side of <- or -> 216 read(c80tmp(:isign-1),"(a)") leftstr 217 ilefttype=0 !close 218 if (index(leftstr,'A')/=0) ilefttype=1 !Alpha 219 if (index(leftstr,'B')/=0) ilefttype=2 !Beta 220 if (ilefttype==0) then 221 read(leftstr,*) orbleft(itmp) 222 else 223 read(leftstr(:len_trim(leftstr)-1),*) orbleft(itmp) 224 if (ilefttype==2) orbleft(itmp)=orbleft(itmp)+nbasis 225 end if 226 !Process right side of <- or -> 227 read(c80tmp(isign+2:),"(a)") rightstr 228 irighttype=0 !close 229 if (index(rightstr,'A')/=0) irighttype=1 !Alpha 230 if (index(rightstr,'B')/=0) irighttype=2 !Beta 231 if (irighttype==0) then 232 read(rightstr,*) orbright(itmp),exccoeff(itmp) 233 else 234 do isplit=1,80 235 if (rightstr(isplit:isplit)=='A'.or.rightstr(isplit:isplit)=='B') exit 236 end do 237 read(rightstr(:isplit-1),*) orbright(itmp) 238 read(rightstr(isplit+1:),*) exccoeff(itmp) 239 if (irighttype==2) orbright(itmp)=orbright(itmp)+nbasis 240 end if 241 end do 242end if 243 244!Test sum of square of the coefficients 245sumsqrexc=0 246sumsqrdeexc=0 247do itmp=1,nexcitorb 248 if (excdir(itmp)==1) sumsqrexc=sumsqrexc+exccoeff(itmp)**2 249 if (excdir(itmp)==2) sumsqrdeexc=sumsqrdeexc-exccoeff(itmp)**2 250end do 251sumsqrall=sumsqrexc+sumsqrdeexc 252write(*,"(' The sum of square of excitation coefficients:',f10.6)") sumsqrexc 253write(*,"(' The negative of the sum of square of de-excitation coefficients:',f10.6)") sumsqrdeexc 254write(*,"(' The sum of above two values',f10.6)") sumsqrall 255close(10) 256exccoeffbackup=exccoeff 257 258 2591 do while(.true.) 260 if (itype==1) then 261 write(*,*) 262 if (idomag==0) write(*,*) "-1 Toggle calculating transit. magnetic dip. density in option 1, current: No" 263 if (idomag==1) write(*,*) "-1 Toggle calculating transit. magnetic dip. density in option 1, current: Yes" 264 write(*,*) "0 Return" 265 write(*,"(a)") " 1 Visualize and analyze hole, electron and transition density and so on" 266 write(*,*) "2 Show contribution of MO pairs to transition dipole moment" 267 write(*,*) "3 Show contribution of each MO to hole and electron distribution" 268 write(*,*) "4 Generate and export transition density matrix" 269 write(*,*) "5 Decompose transition dipole moment to basis function and atom contributions" 270 write(*,*) "6 Calculate Mulliken atomic transition charges" 271 !Note: 7 is not utilized 272 write(*,*) "10 Modify or check excitation coefficients" 273 read(*,*) isel 274 else if (itype==2) then 275 call calcdelta_r(nexcitorb,orbleft,orbright,excdir,exccoeff) 276 return 277 else if (itype==3) then 278 call NTO(nexcitorb,orbleft,orbright,excdir,exccoeff) 279 return 280 else if (itype==4) then 281 call excfragCT(nexcitorb,orbleft,orbright,excdir,exccoeff) 282 return 283 end if 284 285 if (isel==-1) then 286 if (idomag==0) then 287 idomag=1 288 else if (idomag==1) then 289 idomag=0 290 end if 291 else if (isel==0) then 292 return 293 else if (isel==1) then 294 exit 295 !Show contribution of MO pairs to transition dipole moment 296 else if (isel==2) then 297 write(*,"(a)") " Input the threshold for printing, e.g. 0.01 means the orbital pairs having contribution & 298 to any component of transition dipole moment >= 0.01 will be shown" 299 read(*,*) printthres 300 write(*,"(a)") " Input the threshold for calculating contribution, e.g. 0.005 means the configurations with & 301 absolute value of coefficient smaller than 0.005 will be ignored. The smaller the value, the higher the computational cost" 302 read(*,*) calcthres 303 write(*,*) "Generating dipole moment integral matrix..." 304 nsize=nprims*(nprims+1)/2 305 allocate(GTFdipint(3,nsize)) 306 call genGTFDmat(GTFdipint,nsize) 307 allocate(dipcontri(3,nexcitorb)) 308 dipcontri=0 309 fac=1 310 if (wfntype==0.or.wfntype==3) fac=2 311 write(*,*) "Calculating contribution to transition dipole moment, please wait..." 312nthreads=getNThreads() 313!$OMP PARALLEL DO SHARED(dipcontri) PRIVATE(iGTF,jGTF,ides,iexcitorb,imo,jmo) schedule(dynamic) NUM_THREADS(nthreads) 314 do iexcitorb=1,nexcitorb 315 if (abs(exccoeff(iexcitorb))<calcthres) cycle 316! write(*,*) iexcitorb,nexcitorb 317 imo=orbleft(iexcitorb) 318 jmo=orbright(iexcitorb) 319 do iGTF=1,nprims 320 do jGTF=1,nprims 321 if (iGTF>=jGTF) then 322 ides=iGTF*(iGTF-1)/2+jGTF 323 else 324 ides=jGTF*(jGTF-1)/2+iGTF 325 end if 326 dipcontri(:,iexcitorb)=dipcontri(:,iexcitorb)+co(imo,iGTF)*co(jmo,jGTF)*GTFdipint(:,ides) 327 end do 328 end do 329 dipcontri(:,iexcitorb)=dipcontri(:,iexcitorb)*exccoeff(iexcitorb)*fac 330 end do 331!$OMP END PARALLEL DO 332 deallocate(GTFdipint) 333 xdipsum=sum(dipcontri(1,:)) 334 ydipsum=sum(dipcontri(2,:)) 335 zdipsum=sum(dipcontri(3,:)) 336 ishownpair=0 337 write(*,*) " #Pair Coefficient Transition dipole X/Y/Z (au)" 338 do iexcitorb=1,nexcitorb 339 if ( any(abs(dipcontri(:,iexcitorb))>printthres) ) then 340 ishownpair=ishownpair+1 341 imo=orbleft(iexcitorb) 342 jmo=orbright(iexcitorb) 343 strtmp1=" ->" 344 if (excdir(iexcitorb)==2) strtmp1=" <-" 345 if (wfntype==0.or.wfntype==3) then 346 write(*,"(i8,i7,a,i7,f12.6,3f11.6)") iexcitorb,imo,trim(strtmp1),jmo,exccoeff(iexcitorb),dipcontri(:,iexcitorb) 347 else 348 strtmp2="A" 349 strtmp3="A" 350 if (imo>nbasis) then 351 imo=imo-nbasis 352 strtmp2="B" 353 end if 354 if (jmo>nbasis) then 355 jmo=jmo-nbasis 356 strtmp3="B" 357 end if 358 write(*,"(i8,i6,a,a,i6,a,f12.6,3f11.6)") iexcitorb,imo,trim(strtmp2),trim(strtmp1),jmo,trim(strtmp3),exccoeff(iexcitorb),& 359 dipcontri(:,iexcitorb) 360 end if 361 end if 362 end do 363 if (printthres==0) then 364 write(*,"(' Sum: ',3f11.6)") xdipsum,ydipsum,zdipsum 365 else 366 write(*,"(' Sum (including the ones not shown): ',3f11.6)") xdipsum,ydipsum,zdipsum 367! write(*,"(i8,' orbital pairs are shown above')") ishownpair 368 end if 369 oscillstr=2D0/3D0*excenergy/au2eV*(xdipsum**2+ydipsum**2+zdipsum**2) 370 write(*,"(' Norm of total transition dipole moment:',f11.6,' a.u.')") dsqrt(xdipsum**2+ydipsum**2+zdipsum**2) 371 write(*,"(' Oscillator strength:',f12.7)") oscillstr 372 if ((naelec==nbelec).and.iexcmulti==3) write(*,"(a)") " Note: Since the spin multiplicity between the ground state and excited state is different, & 373 the transition dipole moment and thus oscillator strength should be exactly zero in principle" 374 375 write(*,*) 376 write(*,*) "If output all orbital pairs to transdip.txt in current folder? (y/n)" 377 read(*,*) c80tmp 378 if (c80tmp(1:1)=='y'.or.c80tmp(1:1)=='Y') then 379 open(10,file="transdip.txt",status="replace") 380 write(10,*) " #Pair Coefficient Transition dipole X/Y/Z (au)" 381 do iexcitorb=1,nexcitorb 382 imo=orbleft(iexcitorb) 383 jmo=orbright(iexcitorb) 384 strtmp1=" ->" 385 if (excdir(iexcitorb)==2) strtmp1=" <-" 386 if (wfntype==0.or.wfntype==3) then 387 write(10,"(i8,i7,a,i7,f12.6,3f11.6)") iexcitorb,imo,trim(strtmp1),jmo,exccoeff(iexcitorb),dipcontri(:,iexcitorb) 388 else 389 strtmp2="A" 390 strtmp3="A" 391 if (imo>nbasis) then 392 imo=imo-nbasis 393 strtmp2="B" 394 end if 395 if (jmo>nbasis) then 396 jmo=jmo-nbasis 397 strtmp3="B" 398 end if 399 write(10,"(i8,i6,a,a,i6,a,f12.6,3f11.6)") iexcitorb,imo,trim(strtmp2),trim(strtmp1),jmo,trim(strtmp3),exccoeff(iexcitorb),dipcontri(:,iexcitorb) 400 end if 401 end do 402 write(10,"(' Sum: ',3f11.6)") xdipsum,ydipsum,zdipsum 403 write(10,"(' Norm of transition dipole moment:',f12.7,' a.u.')") dsqrt(xdipsum**2+ydipsum**2+zdipsum**2) 404 write(10,"(' Oscillator strength:',f12.7)") oscillstr 405 close(10) 406 write(*,*) "Done! Output finished" 407 end if 408 deallocate(dipcontri) 409 410 !Show contribution of each MO to hole and electron distribution 411 else if (isel==3) then 412 write(*,*) "Input printing criterion" 413 write(*,*) "e.g. 0.005 means only the MOs having contribution >=0.005 will be printed" 414 read(*,*) printcrit 415 allocate(tmparr1(nmo),tmparr2(nmo)) 416 tmparr1=0 !Record contribution of each MO to hole 417 tmparr2=0 !Record contribution of each MO to electron 418 do iexcitorb=1,nexcitorb !Cycle each excitation pair 419 imo=orbleft(iexcitorb) 420 jmo=orbright(iexcitorb) 421 excwei=exccoeff(iexcitorb)**2 422 if (excdir(iexcitorb)==1) then ! -> 423 tmparr1(imo)=tmparr1(imo)+excwei 424 tmparr2(jmo)=tmparr2(jmo)+excwei 425 else ! <- 426 tmparr1(imo)=tmparr1(imo)-excwei 427 tmparr2(jmo)=tmparr2(jmo)-excwei 428 end if 429 end do 430 if (wfntype==0.or.wfntype==3) then !Ground state is close-shell 431 tmparr1=tmparr1*2 432 tmparr2=tmparr2*2 433 do imo=1,nmo 434 if (tmparr1(imo)>=printcrit.or.tmparr2(imo)>=printcrit) & 435 write(*,"(' MO',i8,', Occ:',f10.5,' Hole:',f9.5,' Electron:',f9.5)") imo,MOocc(imo),tmparr1(imo),tmparr2(imo) 436 end do 437 else !Ground state is open-shell 438 do imo=1,nmo 439 if (tmparr1(imo)>=printcrit.or.tmparr2(imo)>=printcrit) then 440 if (imo<=nbasis) then 441 write(*,"(' MO',i7,'A, Occ:',f10.5,' Hole:',f9.5,' Electron:',f9.5)") imo,MOocc(imo),tmparr1(imo),tmparr2(imo) 442 else 443 write(*,"(' MO',i7,'B, Occ:',f10.5,' Hole:',f9.5,' Electron:',f9.5)") imo-nbasis,MOocc(imo),tmparr1(imo),tmparr2(imo) 444 end if 445 end if 446 end do 447 end if 448 write(*,"(' Sum of hole:',f9.5,' Sum of electron:',f9.5)") sum(tmparr1),sum(tmparr2) 449 deallocate(tmparr1,tmparr2) 450 451 !==4: Generating transition density matrix (TDM) 452 !==5: Decompose transition electric/magnetic dipole moment as basis function and atom contributions 453 !==6: Calculate Mulliken atomic transition charges 454 else if (isel==4.or.isel==5.or.isel==6) then 455 if (isel==5) then 456 write(*,*) "Decompose which type of transition dipole moment?" 457 write(*,*) "1: Electric 2: Magnetic" 458 read(*,*) idecomptype 459 end if 460 !There are two ways to construct TDM for TD method, see eqs. 22~24 in JCP,66,3460 for detail. 461 !The first one is correct for transition electric dipole moment, excitation and de-excitation are not distinguished 462 if (isel==4.or.(isel==5.and.idecomptype==1).or.isel==6) iTDMtype=1 463 !The second one is correct for transition velocity/magnetic dipole moment, excitation and de-excitation are considered individually 464 if (isel==5.and.idecomptype==2) iTDMtype=2 465 write(*,*) "Generating transition density matrix..." 466 if (.not.allocated(tdmata)) allocate(tdmata(nbasis,nbasis)) 467 if ((wfntype==1.or.wfntype==4).and.(.not.allocated(tdmatb))) allocate(tdmatb(nbasis,nbasis)) 468 iprog=0 469nthreads=getNThreads() 470!$OMP parallel do shared(tdmata,tdmatb,iprog) private(ibas,jbas,iexcitorb,imo,jmo,tmpval,tmpval2) num_threads(nthreads) SCHEDULE(DYNAMIC) 471 do ibas=1,nbasis 472 do jbas=1,nbasis 473 tmpval=0 474 tmpval2=0 475 if (iTDMtype==1) then 476 do iexcitorb=1,nexcitorb 477 imo=orbleft(iexcitorb) 478 jmo=orbright(iexcitorb) 479 if (MOtype(imo)==0.or.MOtype(imo)==1) then !Close-shell or alpha part of open-shell 480 tmpval=tmpval+exccoeff(iexcitorb)*cobasa(ibas,imo)*cobasa(jbas,jmo) 481 else !Beta part of open-shell 482 tmpval2=tmpval2+exccoeff(iexcitorb)*cobasb(ibas,imo-nbasis)*cobasb(jbas,jmo-nbasis) 483 end if 484 end do 485 else if (iTDMtype==2) then 486 do iexcitorb=1,nexcitorb 487 imo=orbleft(iexcitorb) 488 jmo=orbright(iexcitorb) 489 if (MOtype(imo)==0.or.MOtype(imo)==1) then !Close-shell or alpha part of open-shell 490 if (excdir(iexcitorb)==1) tmpval=tmpval+exccoeff(iexcitorb)*cobasa(ibas,imo)*cobasa(jbas,jmo) 491 if (excdir(iexcitorb)==2) tmpval=tmpval-exccoeff(iexcitorb)*cobasa(ibas,imo)*cobasa(jbas,jmo) 492 else !Beta part of open-shell 493 if (excdir(iexcitorb)==1) tmpval2=tmpval2+exccoeff(iexcitorb)*cobasb(ibas,imo-nbasis)*cobasb(jbas,jmo-nbasis) 494 if (excdir(iexcitorb)==2) tmpval2=tmpval2-exccoeff(iexcitorb)*cobasb(ibas,imo-nbasis)*cobasb(jbas,jmo-nbasis) 495 end if 496 end do 497 end if 498 tdmata(ibas,jbas)=tmpval 499 if (wfntype==1.or.wfntype==4) tdmatb(ibas,jbas)=tmpval2 500 end do 501!$OMP CRITICAL 502 iprog=iprog+1 503 if (nbasis>800) write(*,"(' Progress:',i6,' /',i6)") iprog,nbasis 504!$OMP END CRITICAL 505 end do 506!$OMP END PARALLEL do 507 508 if (wfntype==0.or.wfntype==3) tdmata=tdmata*2 !Close-shell, should double the TDM 509 510 !! Below codes are used to check transition properties based on transition density matrix and corresponding integral matrix 511 !BEWARE THAT CARTESIAN BASIS FUNCTIONS MUST BE USED, SO DON'T FORGET 6D 10F KEYWORDS IN DUE CASES! (However, pure type is OK if you have 512 !set igenDbas/igenMagbas in settings.ini, because the Cartesian integral matrix generated when loading has already been converted to pure case) 513 !Check transition eletric dipole moment. The result is correct when iTDMtype==1, see above 514! if (.not.allocated(Dbas)) then 515! allocate(Dbas(3,nbasis,nbasis)) 516! call genDbas 517! end if 518! Teledipx=sum(tdmata*Dbas(1,:,:)) 519! Teledipy=sum(tdmata*Dbas(2,:,:)) 520! Teledipz=sum(tdmata*Dbas(3,:,:)) 521! write(*,"(' Transition electric dipole moment:',3f12.6)") Teledipx,Teledipy,Teledipz 522! !Check transition velocity dipole moment. The result is correct when iTDMtype==2, see above 523! if (.not.allocated(Velbas)) then 524! allocate(Velbas(3,nbasis,nbasis)) 525! call genvelbas 526! end if 527! Tvdipx=sum(tdmata*Velbas(1,:,:)) 528! Tvdipy=sum(tdmata*Velbas(2,:,:)) 529! Tvdipz=sum(tdmata*Velbas(3,:,:)) 530! write(*,"(' Transition velocity dipole moment:',3f12.6)") Tvdipx,Tvdipy,Tvdipz 531! !Check transition magnetic dipole moment. The result is correct when iTDMtype==2, see above 532! if (.not.allocated(Magbas)) then 533! allocate(Magbas(3,nbasis,nbasis)) 534! call genMagbas 535! end if 536! call showmatgau(tdmata,"tdmata matrix",1) 537! Tmagdipx=sum(tdmata*Magbas(1,:,:)) 538! Tmagdipy=sum(tdmata*Magbas(2,:,:)) 539! Tmagdipz=sum(tdmata*Magbas(3,:,:)) 540! write(*,"(' Transition magnetic dipole moment:',3f12.6)") Tmagdipx,Tmagdipy,Tmagdipz 541! !Below formula is absolutely correct, however the printed result is not correct here because 542! !we cannot obtain correct transition electric and magnetic dipole moments based on the same TDM. 543! !If we directly use the value outputted by Gaussian, you will find the result is in line with the rotatory strength printed by Gaussian 544! !Four points: 545! !1) The negative sign: Because we ignored the negative sign when evaluating magnetic integrals 546! !2) Diveded by two: Necessary by definition 547! !3) 2.54174619D-018: convert electric dipole moment from a.u. to cgs, see gabedit build-in converter 548! !4) 1.85480184D-020: convert magnetic dipole moment from a.u. to cgs, see gabedit build-in converter 549! Rlen=-(Teledipx*Tmagdipx+Teledipy*Tmagdipy+Teledipz*Tmagdipz)/2D0*2.54174619D-018*1.85480184D-020 *1D40 550! write(*,"(' Rotatory strength in length representation:',f14.8,' 10^-40 cgs')") Rlen 551! cycle 552 553 if (isel==4) then !Export transition density matrix 554 write(*,*) "If symmetrizing the transition density matrix? (y/n)" 555 read(*,*) c80tmp 556 if (c80tmp(1:1)=='y'.or.c80tmp(1:1)=='Y') then 557 do ibas=1,nbasis 558 do jbas=ibas,nbasis 559 tdmata(ibas,jbas)=(tdmata(ibas,jbas)+tdmata(jbas,ibas))/dsqrt(2D0) 560 tdmata(jbas,ibas)=tdmata(ibas,jbas) 561 if (wfntype==1.or.wfntype==4) then 562 tdmatb(ibas,jbas)=(tdmatb(ibas,jbas)+tdmatb(jbas,ibas))/dsqrt(2D0) 563 tdmatb(jbas,ibas)=tdmatb(ibas,jbas) 564 end if 565 end do 566 end do 567 end if 568 !Diagonalize transition density matrix, only for test purpose 569! call diagsymat(tdmata,eigvecmat,eigvalarr,istat) 570! call diaggemat(tdmata,eigvecmat,eigvalarr,istat) 571! write(*,"(f12.6)") eigvalarr(:) 572! write(*,*) sum(eigvalarr),sum(eigvalarr,eigvalarr>0),sum(eigvalarr,eigvalarr<0) 573 if (wfntype==0.or.wfntype==3) then 574 open(10,file="tdmat.txt",status="replace") 575 call showmatgau(tdmata,"Transition density matrix",0,"f14.8",10) 576 close(10) 577 write(*,"(a)") " Done! The transition density matrix has been outputted to tdmat.txt in current folder" 578 write(*,"(a)") " Hint: You can use function 2 of main function 18 to plot the transition density matrix by using this file as input file" 579 else 580 open(10,file="tdmata.txt",status="replace") 581 call showmatgau(tdmata,"Alpha transition density matrix",0,"f14.8",10) 582 close(10) 583 open(10,file="tdmatb.txt",status="replace") 584 call showmatgau(tdmatb,"Beta transition density matrix",0,"f14.8",10) 585 close(10) 586 write(*,"(a)") " Done! The alpha and beta transition density matrix has been outputted to tdmata.txt and tdmatb.txt in current folder, respectively" 587 write(*,"(a)") " Hint: You can use function 2 of main function 18 to plot the transition density matrix by using these files as input file" 588 end if 589 590 else if (isel==5) then !Decompose transition electric/magnetic dipole moment as basis function and atom contributions 591 if (idecomptype==1.and.(.not.allocated(Dbas))) then 592 write(*,"(a)") " ERROR: Electric dipole moment integral matrix is not presented. You should set parameter ""igenDbas"" in settings.ini to 1, & 593 so that the matrix will be generated when loading wavefunction file" 594 else if (idecomptype==2.and.(.not.allocated(Magbas))) then 595 write(*,"(a)") " ERROR: Magnetic dipole moment integral matrix is not presented. You should set parameter ""igenMagbas"" in settings.ini to 1, & 596 so that the matrix will be generated when loading wavefunction file" 597 else 598 allocate(trdipmatbas(3,nbasis,nbasis),bastrdip(3,nbasis),trdipmatatm(3,ncenter,ncenter)) 599 if (idecomptype==1) then 600 if (wfntype==1.or.wfntype==4) then 601 write(*,*) "Analyze alpha or beta part? 1=Alpha 2=Beta" 602 read(*,*) iseltmp 603 if (iseltmp==1) then 604 trdipmatbas(1,:,:)=tdmata(:,:)*Dbas(1,:,:) 605 trdipmatbas(2,:,:)=tdmata(:,:)*Dbas(2,:,:) 606 trdipmatbas(3,:,:)=tdmata(:,:)*Dbas(3,:,:) 607 else 608 trdipmatbas(1,:,:)=tdmatb(:,:)*Dbas(1,:,:) 609 trdipmatbas(2,:,:)=tdmatb(:,:)*Dbas(2,:,:) 610 trdipmatbas(3,:,:)=tdmatb(:,:)*Dbas(3,:,:) 611 end if 612 else !Close-shell 613 trdipmatbas(1,:,:)=tdmata(:,:)*Dbas(1,:,:) 614 trdipmatbas(2,:,:)=tdmata(:,:)*Dbas(2,:,:) 615 trdipmatbas(3,:,:)=tdmata(:,:)*Dbas(3,:,:) 616 end if 617 else if (idecomptype==2) then 618 if (wfntype==1.or.wfntype==4) then 619 write(*,*) "Analyze alpha or beta part? 1=Alpha 2=Beta" 620 read(*,*) iseltmp 621 if (iseltmp==1) then 622 trdipmatbas(1,:,:)=tdmata(:,:)*Magbas(1,:,:) 623 trdipmatbas(2,:,:)=tdmata(:,:)*Magbas(2,:,:) 624 trdipmatbas(3,:,:)=tdmata(:,:)*Magbas(3,:,:) 625 else 626 trdipmatbas(1,:,:)=tdmatb(:,:)*Magbas(1,:,:) 627 trdipmatbas(2,:,:)=tdmatb(:,:)*Magbas(2,:,:) 628 trdipmatbas(3,:,:)=tdmatb(:,:)*Magbas(3,:,:) 629 end if 630 else !Close-shell 631 trdipmatbas(1,:,:)=tdmata(:,:)*Magbas(1,:,:) 632 trdipmatbas(2,:,:)=tdmata(:,:)*Magbas(2,:,:) 633 trdipmatbas(3,:,:)=tdmata(:,:)*Magbas(3,:,:) 634 end if 635 end if 636 open(10,file="trdipcontri.txt",status="replace") 637 ides=10 638 write(ides,"(a)") " Contribution from basis functions to transition dipole moment (X,Y,Z, in a.u.) derived by Mulliken partition:" 639 bastrdip=0 640 do ibas=1,nbasis !Note that trdipmatbas is not strictly a symmetry matrix, so we get average value for off-diagonal elements 641 do jbas=1,nbasis 642 if (ibas==jbas) then 643 bastrdip(:,ibas)=bastrdip(:,ibas)+trdipmatbas(:,ibas,jbas) 644 else 645 bastrdip(:,ibas)=bastrdip(:,ibas)+(trdipmatbas(:,ibas,jbas)+trdipmatbas(:,jbas,ibas))/2 646 end if 647 end do 648 write(ides,"(i5,'# Shell:',i5,' Cen:',i5,'(',a2,') Type:',a,3f11.5)") & 649 ibas,basshell(ibas),bascen(ibas),a(bascen(ibas))%name,GTFtype2name(bastype(ibas)),bastrdip(:,ibas) 650 end do 651 write(ides,*) 652 write(ides,"(a)") " Contribution from atoms to transition dipole moment (X,Y,Z, in a.u.) derived by Mulliken partition:" 653 do iatm=1,ncenter 654 write(ides,"(i5,'(',a2,') :',3f12.5)") iatm,a(iatm)%name,& 655 sum(bastrdip(1,basstart(iatm):basend(iatm))),sum(bastrdip(2,basstart(iatm):basend(iatm))),sum(bastrdip(3,basstart(iatm):basend(iatm))) 656 end do 657 write(ides,*) 658 write(ides,"(a,3f12.6,a)") " Transition dipole moment in X/Y/Z",sum(bastrdip(1,:)),sum(bastrdip(2,:)),sum(bastrdip(3,:))," a.u." 659 close(10) 660 write(*,*) "Done! The result has been outputted to trdipcontri.txt in current folder" 661 662 write(*,*) 663 write(*,"(a)") " Would you also like to output atom-atom contribution matrix to AAtrdip.txt in current folder? (y/n)" 664 read(*,*) selectyn 665 if (selectyn=='y') then 666 do iatm=1,ncenter 667 do jatm=1,ncenter 668 trdipmatatm(1,iatm,jatm)=sum( trdipmatbas(1,basstart(iatm):basend(iatm),basstart(jatm):basend(jatm)) ) 669 trdipmatatm(2,iatm,jatm)=sum( trdipmatbas(2,basstart(iatm):basend(iatm),basstart(jatm):basend(jatm)) ) 670 trdipmatatm(3,iatm,jatm)=sum( trdipmatbas(3,basstart(iatm):basend(iatm),basstart(jatm):basend(jatm)) ) 671 !Use the first slot to record sum of square 672 trdipmatatm(1,iatm,jatm)=trdipmatatm(1,iatm,jatm)**2+trdipmatatm(2,iatm,jatm)**2+trdipmatatm(3,iatm,jatm)**2 673 end do 674 end do 675 open(10,file="AAtrdip.txt",status="replace") 676 call showmatgau(trdipmatatm(1,:,:),"Atom-Atom contribution matrix",0,"f14.8",10) 677 close(10) 678 write(*,"(a)") " Done! The matrix has been outputted to AAtrdip.txt in current folder" 679 write(*,"(a)") " This file can be plotted as colored matrix via subfunction 2 of main function 18" 680 end if 681 deallocate(trdipmatbas,bastrdip,trdipmatatm) 682 end if 683 684 else if (isel==6) then !Calculate Mulliken atomic transition charges 685 allocate(bastrpopa(nbasis)) 686 bastrpopa=0 687 if (wfntype==1.or.wfntype==4) then 688 allocate(bastrpopb(nbasis)) 689 bastrpopb=0 690 open(10,file="atmtrchga.chg",status="replace") 691 open(11,file="atmtrchgb.chg",status="replace") 692 else 693 open(10,file="atmtrchg.chg",status="replace") 694 end if 695 do ibas=1,nbasis !Note that tdmat is not strictly a symmetry matrix, so we get average value for off-diagonal elements 696 do jbas=1,nbasis 697 if (ibas==jbas) then 698 bastrpopa(ibas)=bastrpopa(ibas)+tdmata(ibas,jbas)*Sbas(ibas,jbas) 699 if (wfntype==1.or.wfntype==4) bastrpopb(ibas)=bastrpopb(ibas)+tdmatb(ibas,jbas)*Sbas(ibas,jbas) 700 else 701 bastrpopa(ibas)=bastrpopa(ibas)+(tdmata(ibas,jbas)+tdmata(jbas,ibas))*Sbas(ibas,jbas)/2 702 if (wfntype==1.or.wfntype==4) bastrpopb(ibas)=bastrpopb(ibas)+(tdmatb(ibas,jbas)+tdmatb(jbas,ibas))*Sbas(ibas,jbas)/2 703 end if 704 end do 705 end do 706 do iatm=1,ncenter 707 write(10,"(1x,a,4f12.6)") a(iatm)%name,a(iatm)%x*b2a,a(iatm)%y*b2a,a(iatm)%z*b2a,-sum(bastrpopa(basstart(iatm):basend(iatm))) 708 if (wfntype==1.or.wfntype==4) write(11,"(1x,a,4f12.6)") a(iatm)%name,a(iatm)%x*b2a,a(iatm)%y*b2a,a(iatm)%z*b2a,-sum(bastrpopb(basstart(iatm):basend(iatm))) 709 end do 710 close(10) 711 if (wfntype==1.or.wfntype==4) then 712 close(11) 713 write(*,"(a)") " Done! Mulliken atomic transition charges for alpha and beta electrons have been outputted to atmtrchga.chg and atmtrchgb.chg in current folder, respectively" 714 else 715 write(*,"(a)") " Done! Mulliken atomic transition charges have been outputted to atmtrchg.chg in current folder" 716 end if 717 deallocate(bastrpopa) 718 if (wfntype==1.or.wfntype==4) deallocate(bastrpopb) 719 end if 720 deallocate(tdmata) 721 if (wfntype==1.or.wfntype==4) deallocate(tdmatb) 722 723 !Modify or check excitation coefficients 724 else if (isel==10) then 725 write(*,"(' Note: There are',i7,' orbital pairs')") nexcitorb 726 do while(.true.) 727 sumsqrexc=0 728 sumsqrdeexc=0 729 do itmp=1,nexcitorb 730 if (excdir(itmp)==1) sumsqrexc=sumsqrexc+exccoeff(itmp)**2 731 if (excdir(itmp)==2) sumsqrdeexc=sumsqrdeexc-exccoeff(itmp)**2 732 end do 733 write(*,*) 734 write(*,"(' Sum of current coeff.^2 of excitation and de-excitation:',2f11.6)") sumsqrexc,sumsqrdeexc 735 nzerocoeff=count(exccoeff==0) 736 if (nzerocoeff>0) write(*,"(' The number of pairs having zero coefficient:',i7)") nzerocoeff 737! if (wfntype==0.or.wfntype==3) then 738! write(*,"(' The MO index ranges from',i7,' to',i7,', the first',i7,' are occupied')") 1,nmo,nint(naelec) 739! else 740! write(*,"(' The alpha MO index ranges from',i7,' to',i7,', the first',i7,' are occupied')") 1,nbasis,nint(naelec) 741! write(*,"(' The beta MO index ranges from ',i7,' to',i7,', the first',i7,' are occupied')") 1,nbasis,nint(nbelec) 742! end if 743 write(*,*) 744 write(*,*) "-2 Print coefficient (and contribution) of some orbital pairs" 745 write(*,*) "-1 Reset coefficient of all orbital pairs" 746 write(*,*) "0 Return" 747 write(*,*) "1 Set coefficient of an orbital pair" 748 write(*,*) "2 Set coefficient for specific range of orbital pairs" 749 write(*,*) "3 Clean all coefficients but for a specific orbital pair" 750 read(*,*) isel 751 if (isel==-2) then 752 write(*,*) "Input the threshold for printing, e.g. 0.01" 753 write(*,"(a)") " Note: The orbital pairs whose absolute value of coefficient >= this value will be printed. & 754 If input 0, all pairs will be printed." 755 read(*,*) printthres 756 ntmp=0 757 do iexcitorb=1,nexcitorb 758 if (abs(exccoeff(iexcitorb))<printthres) cycle 759 imo=orbleft(iexcitorb) 760 jmo=orbright(iexcitorb) 761 strtmp1=" ->" 762 contrisign=1 763 if (excdir(iexcitorb)==2) then 764 strtmp1=" <-" 765 contrisign=-1 766 end if 767 if (wfntype==0.or.wfntype==3) then 768 write(*,"(i8,i7,a,i7,' Coeff.:',f10.4,' Contri.:',f10.4,'%')") iexcitorb,imo,trim(strtmp1),& 769 jmo,exccoeff(iexcitorb),contrisign*exccoeff(iexcitorb)**2/0.5D0*100 770 else 771 strtmp2="A" 772 strtmp3="A" 773 if (imo>nbasis) then 774 imo=imo-nbasis 775 strtmp2="B" 776 end if 777 if (jmo>nbasis) then 778 jmo=jmo-nbasis 779 strtmp3="B" 780 end if 781 write(*,"(i8,i6,a,a,i6,a,' Coeff.:',f10.4,' Contri.:',f10.4,'%')") iexcitorb,imo,trim(strtmp2),trim(strtmp1),& 782 jmo,trim(strtmp3),exccoeff(iexcitorb),contrisign*exccoeff(iexcitorb)**2*100 783 end if 784 ntmp=ntmp+1 785 end do 786 write(*,"(' There are',i8,' orbital pairs met the criterion you defined')") ntmp 787 else if (isel==-1) then 788 exccoeff=exccoeffbackup 789 write(*,*) "All coefficents have been reset" 790 else if (isel==0) then 791 exit 792 else if (isel==1.or.isel==3) then 793 if (wfntype==0.or.wfntype==3) then 794 write(*,*) "Input the index of the MOs orginally occupied and unoccupied, e.g. 12,26" 795 read(*,*) iMO,jMO 796 else 797 write(*,*) "Input the index of the MOs orginally occupied and unoccupied" 798 write(*,*) "e.g. 12A,46A or 23B,39B" 799 read(*,*) strtmp1,strtmp2 800 read(strtmp1(1:len_trim(strtmp1)-1),*) iMO 801 read(strtmp2(1:len_trim(strtmp2)-1),*) jMO 802 if (index(strtmp1,'B')/=0) then 803 iMO=iMO+nbasis 804 jMO=jMO+nbasis 805 end if 806 end if 807 idir=1 808 if (any(excdir==2)) then 809 write(*,*) "Which is the type of the orbital pair?" 810 if (isel==1) write(*,*) "1: Excitation (namely ->) 2: De-excitation (namely <-)" 811 if (isel==3) write(*,*) "1: Excitation (namely ->) 2: De-excitation (namely <-) 3: Both" 812 read(*,*) idir 813 end if 814 if (isel==1) then !Set coefficient for an orbital pair 815 do iexcitorb=1,nexcitorb 816 if (orbleft(iexcitorb)==iMO.and.orbright(iexcitorb)==jMO.and.excdir(iexcitorb)==idir) then 817 ipair=iexcitorb 818 exit 819 end if 820 end do 821 if (iexcitorb==nexcitorb+1) then 822 write(*,*) "Warning: Cannot find corresponding orbital pair!" 823 else 824 write(*,"(' Original value is',f11.7,', now input the new value, e.g. 0.0143')") exccoeff(ipair) 825 read(*,*) exccoeff(ipair) 826 end if 827 else if (isel==3) then !Clean all coefficients but conserve one 828 ipair1=0 829 ipair2=0 830 do iexcitorb=1,nexcitorb 831 if (orbleft(iexcitorb)==iMO.and.orbright(iexcitorb)==jMO) then 832 if ( excdir(iexcitorb)==1.and.(idir==1.or.idir==3) ) then 833 ipair1=iexcitorb 834 tmp1=exccoeff(ipair1) 835 else if ( excdir(iexcitorb)==2.and.(idir==2.or.idir==3) ) then 836 ipair2=iexcitorb 837 tmp2=exccoeff(ipair2) 838 end if 839 end if 840 end do 841 if ((idir==1.and.ipair1==0).or.(idir2==2.and.ipair2==0).or.(idir==3.and.ipair1==0.and.ipair2==0)) then 842 write(*,*) "Error: Cannot find corresponding orbital pair!" 843 else 844 exccoeff=0 845 if (idir==1.or.(idir==3.and.ipair1/=0)) exccoeff(ipair1)=tmp1 846 if (idir==2.or.(idir==3.and.ipair2/=0)) exccoeff(ipair2)=tmp2 847 write(*,*) "Done!" 848 end if 849 end if 850 else if (isel==2) then 851 if (wfntype==0.or.wfntype==3) then 852 write(*,*) "Input the index range of the MOs orginally occupied, e.g. 14,17" 853 else 854 write(*,*) "Input the index range of the MOs orginally occupied" 855 write(*,*) "e.g. 12A,14A or 22B,28B" 856 end if 857 write(*,*) "Note: Simply input 0 can select all orginally occupied MOs" 858 read(*,"(a)") c80tmp 859 if (c80tmp(1:1)=='0') then 860 istart=1 861 iend=nmo 862 else 863 if (wfntype==0.or.wfntype==3) then 864 read(c80tmp,*) istart,iend 865 else 866 read(c80tmp,*) strtmp1,strtmp2 867 read(strtmp1(1:len_trim(strtmp1)-1),*) istart 868 read(strtmp2(1:len_trim(strtmp2)-1),*) iend 869 if (index(strtmp1,'B')/=0) then 870 istart=istart+nbasis 871 iend=iend+nbasis 872 end if 873 end if 874 end if 875 876 if (wfntype==0.or.wfntype==3) then 877 write(*,*) "Input the index range of the MOs orginally unoccupied, e.g. 72,93" 878 else 879 write(*,*) "Input the index range of the MOs orginally unoccupied" 880 write(*,*) "e.g. 72A,84A or 62B,98B" 881 end if 882 write(*,*) "Note: Simply input 0 can select all orginally unoccupied MOs" 883 read(*,"(a)") c80tmp 884 if (c80tmp(1:1)=='0') then 885 jstart=1 886 jend=nmo 887 else 888 if (wfntype==0.or.wfntype==3) then 889 read(c80tmp,*) jstart,jend 890 else 891 read(c80tmp,*) strtmp1,strtmp2 892 read(strtmp1(1:len_trim(strtmp1)-1),*) jstart 893 read(strtmp2(1:len_trim(strtmp2)-1),*) jend 894 if (index(strtmp1,'B')/=0) then 895 jstart=jstart+nbasis 896 jend=jend+nbasis 897 end if 898 end if 899 end if 900 idir=1 901 if (any(excdir==2)) then 902 write(*,*) "Which is the type of these orbital pairs?" 903 write(*,*) "1: Excitation (namely ->) 2: De-excitation (namely <-) 3: Any" 904 read(*,*) idir 905 end if 906 write(*,*) "Set the coefficients to which value? e.g. 0.00148" 907 read(*,*) tmpval 908 iset=0 909 do iexcitorb=1,nexcitorb 910 if (orbleft(iexcitorb)>=istart.and.orbleft(iexcitorb)<=iend.and.orbright(iexcitorb)>=jstart.and.orbright(iexcitorb)<=jend) then 911 if (idir==3.or.excdir(iexcitorb)==idir) then 912 exccoeff(iexcitorb)=tmpval 913 iset=iset+1 914 end if 915 end if 916 end do 917 write(*,"(' Coefficient of',i7,' orbital pairs are set to',f12.7)") iset,tmpval 918 end if 919 end do 920 end if 921end do 922 923 924!Below we will calculate grid data 925!Set up grid first 926write(*,*) 927call setgrid(0,igridsel) 928if (allocated(holegrid)) deallocate(holegrid,elegrid,holeeleovlp,transdens,holecross,elecross) 929allocate(holegrid(nx,ny,nz),elegrid(nx,ny,nz),holeeleovlp(nx,ny,nz),transdens(nx,ny,nz),holecross(nx,ny,nz),elecross(nx,ny,nz)) 930holegrid=0D0 931elegrid=0D0 932holecross=0D0 933elecross=0D0 934transdens=0D0 935if (idomag==1) then !Will also calculate transition magnetic dipole moment density 936 if (allocated(magtrdens)) deallocate(magtrdens) 937 allocate(magtrdens(nx,ny,nz,3)) 938 magtrdens=0D0 939end if 940write(*,"(a)") " Note: During the calculation of cross term of electron, the orbital pairs whose absolute value of coefficient <0.01 will be ignored to save time" 941allocate(skippair(nexcitorb)) 942skippair=.false. 943do iexcitorb=1,nexcitorb 944 if (abs(exccoeff(iexcitorb))<0.01D0) skippair(iexcitorb)=.true. 945end do 946write(*,*) "Calculating grid data..." 947call walltime(walltime1) 948CALL CPU_TIME(time_begin) 949ifinish=0 950nthreads=getNThreads() 951!$OMP PARALLEL DO SHARED(ifinish,holegrid,elegrid,transdens,holecross,elecross,magtrdens) & 952!$OMP PRIVATE(i,j,k,tmpx,tmpy,tmpz,orbval,wfnderv,imo,jmo,excwei,iexcitorb,jexcitorb,ileft,jleft,iright,jright,tmpleft,tmpright,idir,jdir,tmpval) & 953!$OMP schedule(dynamic) NUM_THREADS(nthreads) 954do k=1,nz 955 tmpz=orgz+(k-1)*dz 956 do j=1,ny 957 tmpy=orgy+(j-1)*dy 958 do i=1,nx 959 tmpx=orgx+(i-1)*dx 960 if (idomag==1) then !Study transition magnetic dipole moment density requests orbital 1st-derivative 961 call orbderv(2,1,nmo,tmpx,tmpy,tmpz,orbval,wfnderv) 962 else 963 call orbderv(1,1,nmo,tmpx,tmpy,tmpz,orbval) 964 end if 965 !Calculate local term of hole and electron 966 do iexcitorb=1,nexcitorb 967 imo=orbleft(iexcitorb) 968 jmo=orbright(iexcitorb) 969 excwei=exccoeff(iexcitorb)**2 970 if (excdir(iexcitorb)==1) then ! -> 971 holegrid(i,j,k)=holegrid(i,j,k)+excwei*orbval(imo)**2 972 elegrid(i,j,k)=elegrid(i,j,k)+excwei*orbval(jmo)**2 973 else ! <- 974 holegrid(i,j,k)=holegrid(i,j,k)-excwei*orbval(imo)**2 975 elegrid(i,j,k)=elegrid(i,j,k)-excwei*orbval(jmo)**2 976! holegrid(i,j,k)=holegrid(i,j,k)+excwei*orbval(jmo)**2 977! elegrid(i,j,k)=elegrid(i,j,k)+excwei*orbval(imo)**2 978 end if 979 transdens(i,j,k)=transdens(i,j,k)+exccoeff(iexcitorb)*orbval(imo)*orbval(jmo) 980 if (idomag==1) then 981 if (excdir(iexcitorb)==1) then 982 magtrdens(i,j,k,1)=magtrdens(i,j,k,1)+exccoeff(iexcitorb)*orbval(imo)*(tmpy*wfnderv(3,jmo)-tmpz*wfnderv(2,jmo)) 983 magtrdens(i,j,k,2)=magtrdens(i,j,k,2)+exccoeff(iexcitorb)*orbval(imo)*(tmpz*wfnderv(1,jmo)-tmpx*wfnderv(3,jmo)) 984 magtrdens(i,j,k,3)=magtrdens(i,j,k,3)+exccoeff(iexcitorb)*orbval(imo)*(tmpx*wfnderv(2,jmo)-tmpy*wfnderv(1,jmo)) 985 else !The de-excitation has important influence on the transition magnetic dipole moment density, so must be considered explicitly 986 magtrdens(i,j,k,1)=magtrdens(i,j,k,1)-exccoeff(iexcitorb)*orbval(imo)*(tmpy*wfnderv(3,jmo)-tmpz*wfnderv(2,jmo)) 987 magtrdens(i,j,k,2)=magtrdens(i,j,k,2)-exccoeff(iexcitorb)*orbval(imo)*(tmpz*wfnderv(1,jmo)-tmpx*wfnderv(3,jmo)) 988 magtrdens(i,j,k,3)=magtrdens(i,j,k,3)-exccoeff(iexcitorb)*orbval(imo)*(tmpx*wfnderv(2,jmo)-tmpy*wfnderv(1,jmo)) 989 end if 990 end if 991 end do 992 !Calculate cross term of hole and electron 993 do iexcitorb=1,nexcitorb 994 !Below cases are skipped: 995 !i->l,i->l and i->l,i<-l and i<-l,i<-l, since ileft==jleft.and.iright==jright 996 !i->l,j->m, since ileft/=jleft.and.iright/=jright 997 !i->l,i<-m and i<-l,j->l, since excdir(iexcitorb)/=excdir(jexcitorb) 998 !**If i->l,i<-l should be taken into account is unsolved 999 ! Currently only take below cases into account: 1000 ! Cross term of hole (do <i|j>): i->l,j->l substract i<-l,j<-l 1001 ! Cross term of electron (do <l|m>): i->l,i->m substract i<-l,i<-m 1002 if (skippair(iexcitorb) .eqv. .true.) cycle 1003 ileft=orbleft(iexcitorb) 1004 iright=orbright(iexcitorb) 1005 tmpleft=exccoeff(iexcitorb)*orbval(ileft) !Use temporary variable to save the time for locating element 1006 tmpright=exccoeff(iexcitorb)*orbval(iright) 1007 idir=excdir(iexcitorb) 1008 do jexcitorb=1,nexcitorb 1009 if (skippair(jexcitorb) .eqv. .true.) cycle 1010 jleft=orbleft(jexcitorb) 1011 jright=orbright(jexcitorb) 1012 jdir=excdir(jexcitorb) 1013 if (idir/=jdir) cycle 1014 if (ileft==jleft) then !do <l|m> 1015 if (iright==jright) cycle 1016 tmpval=tmpright*exccoeff(jexcitorb)*orbval(jright) !Originally virtual orbital 1017 if (idir==1) then !-> 1018 elecross(i,j,k)=elecross(i,j,k)+tmpval 1019 else !<- 1020 elecross(i,j,k)=elecross(i,j,k)-tmpval 1021 end if 1022 else if (iright==jright) then !do <i|j> 1023 tmpval=tmpleft*exccoeff(jexcitorb)*orbval(jleft) !Originally occupied orbital 1024 if (idir==1) then !-> 1025 holecross(i,j,k)=holecross(i,j,k)+tmpval 1026 else !<- 1027 holecross(i,j,k)=holecross(i,j,k)-tmpval 1028 end if 1029 end if 1030 end do 1031 end do 1032 end do 1033 end do 1034 ifinish=ifinish+1 1035 write(*,"(' Finished:',i5,'/',i5)") ifinish,nz 1036end do 1037!$OMP END PARALLEL DO 1038deallocate(skippair) 1039if (wfntype==0.or.wfntype==3) then !For close-shell wavefunction, the weights are normalized to 0.5 (or say the orbitals are doubly occupied), so correct it 1040 holegrid=holegrid*2 1041 elegrid=elegrid*2 1042 transdens=transdens*2 1043 holecross=holecross*2 1044 elecross=elecross*2 1045 if (idomag==1) magtrdens=magtrdens*2 1046end if 1047!Combine local term and cross term of hole to holegrid, that of electron to elegrid. Then local term will be holegrid-holecross and elegrid-elecross 1048holegrid=holegrid+holecross 1049elegrid=elegrid+elecross 1050! holeeleovlp=holegrid*elegrid*min(abs(holegrid),abs(elegrid)) / max(abs(holegrid),abs(elegrid)) 1051! holeeleovlp=holegrid*elegrid 1052holeeleovlp=min(holegrid,elegrid) 1053 1054CALL CPU_TIME(time_end) 1055call walltime(walltime2) 1056write(*,"(' Totally took up CPU time',f12.2,'s, wall clock time',i10,'s',/)") time_end-time_begin,walltime2-walltime1 1057 1058!Check normalization 1059dvol=dx*dy*dz 1060rnormhole=0 1061rnormele=0 1062rnormovlp=0 1063rinttransdens=0 1064rtransdipx=0 1065rtransdipy=0 1066rtransdipz=0 1067centholex=0 1068centholey=0 1069centholez=0 1070centelex=0 1071centeley=0 1072centelez=0 1073rtransmagx=0 1074rtransmagy=0 1075rtransmagz=0 1076do k=1,nz 1077 tmpz=orgz+(k-1)*dz 1078 do j=1,ny 1079 tmpy=orgy+(j-1)*dy 1080 do i=1,nx 1081 tmpx=orgx+(i-1)*dx 1082 rnormhole=rnormhole+holegrid(i,j,k) 1083 rnormele=rnormele+elegrid(i,j,k) 1084 rnormovlp=rnormovlp+holeeleovlp(i,j,k) 1085 rinttransdens=rinttransdens+transdens(i,j,k) 1086 rtransdipx=rtransdipx-tmpx*transdens(i,j,k) 1087 rtransdipy=rtransdipy-tmpy*transdens(i,j,k) 1088 rtransdipz=rtransdipz-tmpz*transdens(i,j,k) 1089 centholex=centholex+holegrid(i,j,k)*tmpx 1090 centholey=centholey+holegrid(i,j,k)*tmpy 1091 centholez=centholez+holegrid(i,j,k)*tmpz 1092 centelex=centelex+elegrid(i,j,k)*tmpx 1093 centeley=centeley+elegrid(i,j,k)*tmpy 1094 centelez=centelez+elegrid(i,j,k)*tmpz 1095 if (idomag==1) then 1096 rtransmagx=rtransmagx+magtrdens(i,j,k,1) 1097 rtransmagy=rtransmagy+magtrdens(i,j,k,2) 1098 rtransmagz=rtransmagz+magtrdens(i,j,k,3) 1099 end if 1100 end do 1101 end do 1102end do 1103rnormhole=rnormhole*dvol 1104rnormele=rnormele*dvol 1105rnormovlp=rnormovlp*dvol 1106rinttransdens=rinttransdens*dvol 1107rtransdipx=rtransdipx*dvol 1108rtransdipy=rtransdipy*dvol 1109rtransdipz=rtransdipz*dvol 1110centholex=centholex*dvol/rnormhole 1111centholey=centholey*dvol/rnormhole 1112centholez=centholez*dvol/rnormhole 1113centelex=centelex*dvol/rnormele 1114centeley=centeley*dvol/rnormele 1115centelez=centelez*dvol/rnormele 1116rtransmagx=rtransmagx*dvol 1117rtransmagy=rtransmagy*dvol 1118rtransmagz=rtransmagz*dvol 1119write(*,"(' Integral of hole:',f12.6)") rnormhole 1120write(*,"(' Integral of electron:',f12.6)") rnormele 1121write(*,"(' Integral of transition density:',f12.6)") rinttransdens 1122write(*,"(' Transition dipole moment in X/Y/Z:',3f11.6,' a.u.')") rtransdipx,rtransdipy,rtransdipz 1123if (idomag==1) write(*,"(' Transition magnetic dipole moment in X/Y/Z:',3f10.6,' a.u.')") rtransmagx,rtransmagy,rtransmagz 1124write(*,"(' Integral of overlap of hole-electron distribution:',f12.7)") rnormovlp 1125write(*,"(' Centroid of hole in X/Y/Z: ',3f12.6,' Angstrom')") centholex*b2a,centholey*b2a,centholez*b2a 1126write(*,"(' Centroid of electron in X/Y/Z:',3f12.6,' Angstrom')") centelex*b2a,centeley*b2a,centelez*b2a 1127write(*,"(' Distance between centroid of hole and electron in X/Y/Z:')") 1128disx=abs(centelex-centholex) 1129disy=abs(centeley-centholey) 1130disz=abs(centelez-centholez) 1131disnorm=dsqrt(disx**2+disy**2+disz**2) 1132write(*,"(3f12.6,' Angstrom Norm:',f12.6' Angstrom')") disx*b2a,disy*b2a,disz*b2a,disnorm*b2a 1133!Of course, if relaxed density it used, we cannot evaluate the variation of dipole moment here exactly 1134write(*,"(' Variation of dipole moment with respect to ground state in X/Y/Z:')") 1135avgtransval=(rnormele+rnormhole)/2D0 !Ideal value is 1.0 1136dipvarx=-(centelex-centholex)*avgtransval 1137dipvary=-(centeley-centholey)*avgtransval 1138dipvarz=-(centelez-centholez)*avgtransval 1139write(*,"(3f12.6,' a.u. Norm:',f12.6' a.u.')") dipvarx,dipvary,dipvarz,dsqrt(dipvarx**2+dipvary**2+dipvarz**2) 1140 1141sigxele=0D0 1142sigyele=0D0 1143sigzele=0D0 1144sigxhole=0D0 1145sigyhole=0D0 1146sigzhole=0D0 1147do k=1,nz 1148 tmpz=orgz+(k-1)*dz 1149 do j=1,ny 1150 tmpy=orgy+(j-1)*dy 1151 do i=1,nx 1152 tmpx=orgx+(i-1)*dx 1153 sigxele=sigxele+elegrid(i,j,k)*(tmpx-centelex)**2 1154 sigyele=sigyele+elegrid(i,j,k)*(tmpy-centeley)**2 1155 sigzele=sigzele+elegrid(i,j,k)*(tmpz-centelez)**2 1156 sigxhole=sigxhole+holegrid(i,j,k)*(tmpx-centholex)**2 1157 sigyhole=sigyhole+holegrid(i,j,k)*(tmpy-centholey)**2 1158 sigzhole=sigzhole+holegrid(i,j,k)*(tmpz-centholez)**2 1159 end do 1160 end do 1161end do 1162sigxele=dsqrt(sigxele*dvol/rnormele) 1163sigyele=dsqrt(sigyele*dvol/rnormele) 1164sigzele=dsqrt(sigzele*dvol/rnormele) 1165signormele=dsqrt(sigxele**2+sigyele**2+sigzele**2) 1166sigxhole=dsqrt(sigxhole*dvol/rnormhole) 1167sigyhole=dsqrt(sigyhole*dvol/rnormhole) 1168sigzhole=dsqrt(sigzhole*dvol/rnormhole) 1169signormhole=dsqrt(sigxhole**2+sigyhole**2+sigzhole**2) 1170write(*,"(' RMSD of electron in x,y,z:',3f8.3,' Total:',f8.3,' Angstrom')") sigxele*b2a,sigyele*b2a,sigzele*b2a,signormele*b2a 1171write(*,"(' RMSD of hole in x,y,z: ',3f8.3,' Total:',f8.3,' Angstrom')") sigxhole*b2a,sigyhole*b2a,sigzhole*b2a,signormhole*b2a 1172delta_sigCT=dsqrt((sigxele-sigxhole)**2+(sigyele-sigyhole)**2+(sigzele-sigzhole)**2) 1173write(*,"(' Difference between RMSD of hole and electron:',f8.3,' Angstrom')") delta_sigCT*b2a 1174Hx=(sigxele+sigxhole)/2D0 1175Hy=(sigyele+sigyhole)/2D0 1176Hz=(sigzele+sigzhole)/2D0 1177Hnorm=dsqrt(Hx**2+Hy**2+Hz**2) 1178write(*,"(' H index in x,y,z:',3f8.3,' Norm:',f8.3,' Angstrom')") Hx*b2a,Hy*b2a,Hz*b2a,Hnorm*b2a 1179tx=disx-Hx 1180ty=disy-Hy 1181tz=disz-Hz 1182tnorm=dsqrt(tx**2+ty**2+tz**2) 1183write(*,"(' t index in x,y,z:',3f8.3,' Norm:',f8.3,' Angstrom')") tx*b2a,ty*b2a,tz*b2a,tnorm*b2a 1184 1185!---- Calculate ghost state diagnostic index proposed by Adamo (DOI: 10.1002/jcc.24862) 1186ghostp2=1/disnorm !in a.u. 1187!Definition 1 (the original paper definition). The original paper doesn't clearly mention how to deal with de-excitation, so I treat it as usual 1188sumC=sum(exccoeff(1:nexcitorb)) 1189ghostp1=0 1190do itmp=1,nexcitorb 1191 ghostp1=ghostp1+exccoeff(itmp)/sumC*(-MOene(orbleft(itmp))-MOene(orbright(itmp))) 1192end do 1193ghostidx=ghostp1-ghostp2 1194write(*,"(' Ghost-hunter index (defin. 1):',f8.3,' eV, 1st/2nd terms:',2f9.3,' eV')") ghostidx*au2eV,ghostp1*au2eV,ghostp2*au2eV 1195!Definition 2 (defined by Tian Lu, more reasonable) 1196sumCsqr=0 1197do itmp=1,nexcitorb 1198 if (excdir(itmp)==2) cycle !Skip de-excitation 1199 sumCsqr=sumCsqr+exccoeff(itmp)**2 1200end do 1201ghostp1=0 1202do itmp=1,nexcitorb 1203 if (excdir(itmp)==2) cycle !Skip de-excitation 1204 ghostp1=ghostp1+exccoeff(itmp)**2/sumCsqr*(-MOene(orbleft(itmp))-MOene(orbright(itmp))) 1205end do 1206ghostidx=ghostp1-ghostp2 1207write(*,"(' Ghost-hunter index (defin. 2):',f8.3,' eV, 1st/2nd terms:',2f9.3,' eV')") ghostidx*au2eV,ghostp1*au2eV,ghostp2*au2eV 1208write(*,"(' Excitation energy of this state:',f10.3,' eV')") excenergy 1209if (excenergy<ghostidx*au2eV) write(*,*) "Warning: Probably this is a ghost state" 1210 1211if (allocated(Cele)) deallocate(Cele,Chole) 1212allocate(Cele(nx,ny,nz),Chole(nx,ny,nz)) 1213do i=1,nx 1214 do j=1,ny 1215 do k=1,nz 1216 rnowx=orgx+(i-1)*dx 1217 rnowy=orgy+(j-1)*dy 1218 rnowz=orgz+(k-1)*dz 1219 Cele(i,j,k)=exp( -(rnowx-centelex)**2/(2*sigxele**2) -(rnowy-centeley)**2/(2*sigyele**2) -(rnowz-centelez)**2/(2*sigzele**2)) 1220 Chole(i,j,k)=exp( -(rnowx-centholex)**2/(2*sigxhole**2) -(rnowy-centholey)**2/(2*sigyhole**2) -(rnowz-centholez)**2/(2*sigzhole**2)) 1221 end do 1222 end do 1223end do 1224Cele=Cele*rnormele/(sum(Cele)*dvol) 1225Chole=Chole*rnormhole/(sum(Chole)*dvol) 1226 1227do while(.true.) 1228 write(*,*) 1229 write(*,*) "0 Return" 1230 write(*,*) "1 Show isosurface of hole distribution" 1231 write(*,*) "2 Show isosurface of electron distribution" 1232 write(*,*) "3 Show isosurface of hole and electron distribution simultaneously" 1233 write(*,*) "4 Show isosurface of overlap of hole-electron" 1234 write(*,*) "5 Show isosurface of transition density" 1235 write(*,*) "6 Show isosurface of transition dipole moment density" 1236 write(*,*) "7 Show isosurface of charge density difference" 1237 write(*,*) "8 Show isosurface of Cele and Chole functions simultaneously" 1238 if (idomag==1) write(*,*) "9 Show isosurface of transition magnetic dipole moment density" 1239 write(*,*) "10 Output cube file of hole distribution to current folder" 1240 write(*,*) "11 Output cube file of electron distribution to current folder" 1241 write(*,*) "12 Output cube file of overlap of hole-electron to current folder" 1242 write(*,*) "13 Output cube file of transition density to current folder" 1243 write(*,*) "14 Output cube file of transition dipole moment density to current folder" 1244 write(*,*) "15 Output cube file of charge density difference to current folder" 1245 write(*,*) "16 Output cube file of Cele and Chole functions to current folder" 1246 write(*,*) "17 Output cube file of transition magnetic dipole moment density" 1247 write(*,*) "18 Calculate hole-electron Coulomb attractive energy" 1248 read(*,*) isel 1249 if (isel==0) then 1250 goto 1 1251 else if (isel==1.or.isel==2.or.isel==4.or.isel==5.or.isel==6.or.isel==7.or.isel==9) then 1252 if (allocated(cubmat)) deallocate(cubmat) 1253 allocate(cubmat(nx,ny,nz)) 1254 if (isel==1) then 1255 write(*,*) "Select the type of hole distribution to be shown. 1 is commonly used" 1256 write(*,*) "1 Total (local term + cross term)" 1257 write(*,*) "2 Local term only" 1258 write(*,*) "3 Cross term only" 1259 read(*,*) isel2 1260 if (isel2==1) then 1261 cubmat=holegrid 1262 else if (isel2==2) then 1263 cubmat=holegrid-holecross 1264 else if (isel2==3) then 1265 cubmat=holecross 1266 end if 1267 sur_value=0.002D0 1268 else if (isel==2) then 1269 write(*,*) "Select the type of electron distribution to be shown. 1 is commonly used" 1270 write(*,*) "1 Total (local term + cross term)" 1271 write(*,*) "2 Local term only" 1272 write(*,*) "3 Cross term only" 1273 read(*,*) isel2 1274 if (isel2==1) then 1275 cubmat=elegrid 1276 else if (isel2==2) then 1277 cubmat=elegrid-elecross 1278 else if (isel2==3) then 1279 cubmat=elecross 1280 end if 1281 sur_value=0.002D0 1282 else if (isel==4) then 1283 cubmat=holeeleovlp 1284 sur_value=0.002D0 1285 else if (isel==5) then 1286 cubmat=transdens 1287 isosur1style=4 1288 sur_value=0.001D0 1289 else if (isel==6) then 1290 write(*,*) "Select the component of transition dipole moment density" 1291 write(*,*) "1: X component 2: Y component 3: Z component" 1292 read(*,*) ifac 1293 isosur1style=4 1294 sur_value=0.001D0 1295 do k=1,nz 1296 tmpz=orgz+(k-1)*dz 1297 do j=1,ny 1298 tmpy=orgy+(j-1)*dy 1299 do i=1,nx 1300 tmpx=orgx+(i-1)*dx 1301 if (ifac==1) cubmat(i,j,k)=-tmpx*transdens(i,j,k) 1302 if (ifac==2) cubmat(i,j,k)=-tmpy*transdens(i,j,k) 1303 if (ifac==3) cubmat(i,j,k)=-tmpz*transdens(i,j,k) 1304 end do 1305 end do 1306 end do 1307 else if (isel==7) then 1308 cubmat=elegrid-holegrid 1309 sur_value=0.002D0 1310 else if (isel==9) then 1311 write(*,*) "Select the component of transition magnetic dipole moment density" 1312 write(*,*) "1: X component 2: Y component 3: Z component" 1313 read(*,*) ifac 1314 isosur1style=4 1315 sur_value=0.007D0 1316 if (ifac==1) cubmat=magtrdens(:,:,:,1) 1317 if (ifac==2) cubmat=magtrdens(:,:,:,2) 1318 if (ifac==3) cubmat=magtrdens(:,:,:,3) 1319 end if 1320 deallocate(cubmat) 1321 else if (isel==3.or.isel==8) then 1322 if (isel==3) write(*,"(a)") " Note: Blue and green isosurfaces represent hole and electron distributions, respectively" 1323 if (isel==8) write(*,"(a)") " Note: Blue and green isosurfaces represent Chole and Cele functions, respectively" 1324 sur_value=0.002D0 1325 isosursec=1 1326 clrRcub2sameold=clrRcub2same !Backup previous color setting 1327 clrGcub2sameold=clrGcub2same 1328 clrBcub2sameold=clrBcub2same 1329 clrRcub2samemeshptold=clrRcub2samemeshpt 1330 clrGcub2samemeshptold=clrGcub2samemeshpt 1331 clrBcub2samemeshptold=clrBcub2samemeshpt 1332 clrRcub2same=0.3D0 1333 clrGcub2same=0.45D0 1334 clrBcub2same=0.9D0 1335 clrRcub2samemeshpt=0.3D0 1336 clrGcub2samemeshpt=0.45D0 1337 clrBcub2samemeshpt=0.9D0 1338 if (allocated(cubmat)) deallocate(cubmat) 1339 if (allocated(cubmattmp)) deallocate(cubmattmp) 1340 allocate(cubmat(nx,ny,nz),cubmattmp(nx,ny,nz)) 1341 if (isel==3) then 1342 cubmat=elegrid 1343 cubmattmp=holegrid 1344 else if (isel==8) then 1345 cubmat=Cele 1346 cubmattmp=Chole 1347 end if 1348 isosur1style=4 1349 isosur2style=4 1350 deallocate(cubmat,cubmattmp) 1351 clrRcub2same=clrRcub2sameold !Recover previous color setting 1352 clrGcub2same=clrGcub2sameold 1353 clrBcub2same=clrBcub2sameold 1354 clrRcub2samemeshpt=clrRcub2samemeshptold 1355 clrGcub2samemeshpt=clrGcub2samemeshptold 1356 clrBcub2samemeshpt=clrBcub2samemeshptold 1357 else if (isel==10) then 1358 write(*,*) "Select the type of hole distribution to be exported. 1 is commonly used" 1359 write(*,*) "1 Total (local term + cross term)" 1360 write(*,*) "2 Local term only" 1361 write(*,*) "3 Cross term only" 1362 read(*,*) isel2 1363 write(*,*) "Outputting hole distribution to hole.cub in current folder" 1364 open(10,file="hole.cub",status="replace") 1365 if (isel2==1) call outcube(holegrid,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1366 if (isel2==2) call outcube(holegrid-holecross,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1367 if (isel2==3) call outcube(holecross,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1368 close(10) 1369 write(*,*) "Done!" 1370 else if (isel==11) then 1371 write(*,*) "Select the type of electron distribution to be exported. 1 is commonly used" 1372 write(*,*) "1 Total (local term + cross term)" 1373 write(*,*) "2 Local term only" 1374 write(*,*) "3 Cross term only" 1375 read(*,*) isel2 1376 write(*,*) "Outputting electron distribution to electron.cub in current folder" 1377 open(10,file="electron.cub",status="replace") 1378 if (isel2==1) call outcube(elegrid,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1379 if (isel2==2) call outcube(elegrid-elecross,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1380 if (isel2==3) call outcube(elecross,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1381 close(10) 1382 write(*,*) "Done!" 1383 else if (isel==12) then 1384 write(*,"(a)") " Outputting overlap of hole-electron to holeeleovlp.cub in current folder" 1385 open(10,file="holeeleovlp.cub",status="replace") 1386 call outcube(holeeleovlp,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1387 close(10) 1388 write(*,*) "Done!" 1389 else if (isel==13) then 1390 write(*,"(a)") " Outputting transition density to transdens.cub in current folder" 1391 open(10,file="transdens.cub",status="replace") 1392 call outcube(transdens,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1393 close(10) 1394 write(*,*) "Done!" 1395 else if (isel==14) then 1396 write(*,*) "Select the component of transition dipole moment density" 1397 write(*,*) "1: X component 2: Y component 3: Z component" 1398 read(*,*) ifac 1399 if (allocated(cubmat)) deallocate(cubmat) 1400 allocate(cubmat(nx,ny,nz)) 1401 do k=1,nz 1402 tmpz=orgz+(k-1)*dz 1403 do j=1,ny 1404 tmpy=orgy+(j-1)*dy 1405 do i=1,nx 1406 tmpx=orgx+(i-1)*dx 1407 if (ifac==1) cubmat(i,j,k)=-tmpx*transdens(i,j,k) 1408 if (ifac==2) cubmat(i,j,k)=-tmpy*transdens(i,j,k) 1409 if (ifac==3) cubmat(i,j,k)=-tmpz*transdens(i,j,k) 1410 end do 1411 end do 1412 end do 1413 write(*,*) "Outputting to transdipdens.cub in current folder" 1414 open(10,file="transdipdens.cub",status="replace") 1415 call outcube(cubmat,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1416 close(10) 1417 write(*,*) "Done!" 1418 deallocate(cubmat) 1419 else if (isel==15) then 1420 write(*,*) "Outputting charge density difference to CDD.cub in current folder" 1421 open(10,file="CDD.cub",status="replace") 1422 call outcube(elegrid-holegrid,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1423 close(10) 1424 write(*,*) "Done!" 1425 else if (isel==16) then 1426 open(10,file="Cele.cub",status="replace") 1427 call outcube(Cele,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1428 close(10) 1429 write(*,"('Cele function has been outputted to ""Cele.cub"" in current folder')") 1430 open(10,file="Chole.cub",status="replace") 1431 call outcube(Chole,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1432 close(10) 1433 write(*,"('Chole function has been outputted to ""Chole.cub"" in current folder')") 1434 else if (isel==17) then 1435 write(*,*) "Select the component of transition magnetic dipole moment density" 1436 write(*,*) "1: X component 2: Y component 3: Z component" 1437 read(*,*) ifac 1438 write(*,*) "Outputting to magtrdipdens.cub in current folder" 1439 open(10,file="magtrdipdens.cub",status="replace") 1440 if (ifac==1) call outcube(magtrdens(:,:,:,1),nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1441 if (ifac==2) call outcube(magtrdens(:,:,:,2),nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1442 if (ifac==3) call outcube(magtrdens(:,:,:,3),nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 1443 close(10) 1444 write(*,*) "Done!" 1445 else if (isel==18) then 1446 call walltime(iwalltime1) 1447 CALL CPU_TIME(time_begin) 1448 write(*,*) "Calculating, please wait..." 1449 allocate(cubx(nx),cuby(ny),cubz(nz)) 1450 do i=1,nx 1451 cubx(i)=orgx+(i-1)*dx 1452 end do 1453 do i=1,ny 1454 cuby(i)=orgy+(i-1)*dy 1455 end do 1456 do i=1,nz 1457 cubz(i)=orgz+(i-1)*dz 1458 end do 1459 coulene=0 1460 do i=1,nx 1461 do j=1,ny 1462 do k=1,nz 1463 if (Cele(i,j,k)<1D-6) cycle !Typically leads to error at 0.001 magnitude 1464nthreads=getNThreads() 1465!$OMP parallel shared(coulene) private(ii,jj,kk,distx2,disty2,distz2,dist,coulenetmp) num_threads(nthreads) 1466 coulenetmp=0 1467!$OMP do schedule(DYNAMIC) 1468 do ii=1,nx 1469 distx2=(cubx(i)-cubx(ii))**2 1470 do jj=1,ny 1471 disty2=(cuby(j)-cuby(jj))**2 1472 do kk=1,nz 1473 if (i==ii.and.j==jj.and.k==kk) cycle 1474 distz2=(cubz(k)-cubz(kk))**2 1475 dist=dsqrt(distx2+disty2+distz2) 1476 coulenetmp=coulenetmp+Cele(i,j,k)*Chole(ii,jj,kk)/dist 1477 end do 1478 end do 1479 end do 1480!$OMP END DO 1481!$OMP CRITICAL 1482 coulene=coulene+coulenetmp 1483!$OMP END CRITICAL 1484!$OMP END PARALLEL 1485 end do 1486 end do 1487 write(*,"(' Finished:',i4,' /',i4)") i,nx 1488 end do 1489 dvol=dx*dy*dz 1490 coulene=-coulene*dvol*dvol 1491 CALL CPU_TIME(time_end) 1492 call walltime(iwalltime2) 1493 write(*,"(' Calculation took up CPU time',f12.2,'s, wall clock time',i10,'s')") time_end-time_begin,iwalltime2-iwalltime1 1494 write(*,*) 1495 write(*,"(' Coulomb attractive energy:',f12.6,' a.u. (',f12.6,' eV )')") coulene,coulene*au2eV 1496 deallocate(cubx,cuby,cubz) 1497 end if 1498end do 1499end subroutine 1500 1501 1502 1503 1504!!--- Calculate delta_r index, see J. Chem. Theory Comput., 9, 3118 (2013) 1505!excitation and de-excitation cofficients are summed together, according to Eq.9 of the original paper 1506subroutine calcdelta_r(nexcitorb,orbleft,orbright,excdir,exccoeff) 1507use defvar 1508implicit real*8 (a-h,o-z) 1509integer nexcitorb,orbleft(nexcitorb),orbright(nexcitorb),excdir(nexcitorb) 1510real*8 exccoeff(nexcitorb) 1511real*8,allocatable :: exccoefftot(:) !Store the coefficient combined from excitation and de-excitation of the same pair 1512real*8,allocatable :: orbcenx(:),orbceny(:),orbcenz(:) !Store centroid of MOs 1513real*8,allocatable :: GTFdipint(:,:) 1514character strtmp1,strtmp2,selectyn 1515write(*,*) "Calculating, please wait..." 1516allocate(exccoefftot(nexcitorb)) 1517exccoefftot=exccoeff 1518coeffsumsqr=0D0 1519do itmp=1,nexcitorb 1520 if (excdir(itmp)==1) then !->, find corresponding <- pair and sum its coefficient to here 1521 do jtmp=1,nexcitorb 1522 if (excdir(jtmp)==2.and.orbleft(itmp)==orbleft(jtmp).and.orbright(itmp)==orbright(jtmp)) exccoefftot(itmp)=exccoefftot(itmp)+exccoeff(jtmp) 1523 end do 1524 end if 1525 coeffsumsqr=coeffsumsqr+exccoefftot(itmp)**2 1526end do 1527!Calculate dipole moment integral matrix 1528nsize=nprims*(nprims+1)/2 1529allocate(GTFdipint(3,nsize)) 1530call genGTFDmat(GTFdipint,nsize) 1531!Calculate centroid of MOs 1532allocate(orbcenx(nmo),orbceny(nmo),orbcenz(nmo)) 1533orbcenx=0 1534orbceny=0 1535orbcenz=0 1536nthreads=getNThreads() 1537!$OMP PARALLEL DO SHARED(orbcenx,orbceny,orbcenz) PRIVATE(iGTF,jGTF,ides,imo) schedule(dynamic) NUM_THREADS(nthreads) 1538do imo=1,nmo 1539 do iGTF=1,nprims 1540 do jGTF=1,nprims 1541 if (iGTF>=jGTF) then 1542 ides=iGTF*(iGTF-1)/2+jGTF 1543 else 1544 ides=jGTF*(jGTF-1)/2+iGTF 1545 end if 1546 orbcenx(imo)=orbcenx(imo)+co(imo,iGTF)*co(imo,jGTF)*GTFdipint(1,ides) 1547 orbceny(imo)=orbceny(imo)+co(imo,iGTF)*co(imo,jGTF)*GTFdipint(2,ides) 1548 orbcenz(imo)=orbcenz(imo)+co(imo,iGTF)*co(imo,jGTF)*GTFdipint(3,ides) 1549 end do 1550 end do 1551end do 1552!$OMP END PARALLEL DO 1553deallocate(GTFdipint) 1554delta_r=0 1555do itmp=1,nexcitorb 1556 if (excdir(itmp)==2) cycle 1557 imo=orbleft(itmp) 1558 jmo=orbright(itmp) 1559 delta_r=delta_r+exccoefftot(itmp)**2/coeffsumsqr *dsqrt((orbcenx(imo)-orbcenx(jmo))**2+(orbceny(imo)-orbceny(jmo))**2+(orbcenz(imo)-orbcenz(jmo))**2) 1560end do 1561write(*,"(' Delta_r=',f12.6,' Bohr,',f12.6,' Angstrom')") delta_r,delta_r*b2a 1562write(*,*) 1563write(*,*) "If print orbital pair contribution to delta_r? (y/n)" 1564read(*,*) selectyn 1565if (selectyn=='y'.or.selectyn=='Y') then 1566 write(*,*) "Input threshold for printing e.g. 0.05" 1567 write(*,"(a)") " Note: If input -1, then all contributions will be exported to delta_r.txt in curren folder" 1568 read(*,*) printthres 1569 iout=6 1570 if (printthres==-1) then 1571 iout=10 1572 open(10,file="delta_r.txt",status="replace") 1573 end if 1574 write(iout,"(a)") " Note: The transition coefficients shown below have combined both excitation and de-excitation parts" 1575 write(iout,"(' Sum of square of transition coefficient:',f12.6)") coeffsumsqr 1576 write(iout,*) " #Pair Orbitals Coefficient Contribution (Bohr and Angstrom)" 1577 do itmp=1,nexcitorb 1578 if (excdir(itmp)==2) cycle 1579 imo=orbleft(itmp) 1580 jmo=orbright(itmp) 1581 contrival=exccoefftot(itmp)**2/coeffsumsqr *dsqrt((orbcenx(imo)-orbcenx(jmo))**2+(orbceny(imo)-orbceny(jmo))**2+(orbcenz(imo)-orbcenz(jmo))**2) 1582 if (contrival<printthres) cycle 1583 if (wfntype==0.or.wfntype==3) then 1584 write(iout,"(i8,2i7,f16.7,3x,2f16.7)") itmp,imo,jmo,exccoefftot(itmp),contrival,contrival*b2a 1585 else 1586 strtmp1="A" 1587 strtmp2="A" 1588 if (imo>nbasis) then 1589 imo=imo-nbasis 1590 strtmp1="B" 1591 end if 1592 if (jmo>nbasis) then 1593 jmo=jmo-nbasis 1594 strtmp2="B" 1595 end if 1596 write(iout,"(i8,i6,a,i6,a,f16.7,3x,2f16.7)") itmp,imo,strtmp1,jmo,strtmp2,exccoefftot(itmp),contrival,contrival*b2a 1597 end if 1598 end do 1599 if (printthres==-1) then 1600 close(10) 1601 write(*,*) "Done, outputting finished" 1602 end if 1603end if 1604deallocate(exccoefftot,orbcenx,orbceny,orbcenz) 1605write(*,*) 1606end subroutine 1607 1608 1609 1610 1611!!---- Generate NTOs, original paper of NTO: J. Chem. Phys., 118, 4775 (2003) 1612!The NTO eigenvalues for unrestricted wavefunction is 1/2 of the ones outputted by Gaussian, & 1613!which is incorrect (i.e. the sum is 2.0 rather than 1.0), the result must be Gaussian used incorrect factor 1614subroutine NTO(nexcitorb,orbleft,orbright,excdir,exccoeff) 1615use defvar 1616use util 1617implicit real*8 (a-h,o-z) 1618real*8,allocatable :: T_MO(:,:),TT(:,:),NTOvec(:,:),NTOval(:) 1619integer nexcitorb,orbleft(nexcitorb),orbright(nexcitorb),excdir(nexcitorb) 1620real*8 exccoeff(nexcitorb),tmparr(nbasis) 1621character c200tmp*200 1622write(*,*) 1623NTOvalcoeff=2 1624if (allocated(CObasb)) then 1625 NTOvalcoeff=1 1626 write(*,*) "Result of Alpha part:" 1627end if 1628!*** Alpha part or restricted case 1629nocc=nint(naelec) 1630nvir=nbasis-nocc 1631allocate(T_MO(nocc,nvir)) 1632T_MO=0 1633do iexcitorb=1,nexcitorb 1634 iocc=orbleft(iexcitorb) 1635 ivir=orbright(iexcitorb)-nocc 1636 if (iocc>nbasis) cycle !Here we only process alpha part, index of Beta orbitals are higher than nbasis 1637 !Ignoring de-excitation, this treatment makes the result identical to that ouputted by Gaussian 1638 if (excdir(iexcitorb)==1) T_MO(iocc,ivir)=T_MO(iocc,ivir)+exccoeff(iexcitorb) 1639end do 1640!Occupied part 1641allocate(TT(nocc,nocc),NTOvec(nocc,nocc),NTOval(nocc)) 1642TT=matmul(T_MO,transpose(T_MO)) 1643call diagsymat(TT,NTOvec,NTOval,ierror) 1644NTOval=NTOval*NTOvalcoeff 1645MOene(1:nocc)=NTOval !By default, the diagsymat gives result from low to high 1646CObasa(:,1:nocc)=matmul(CObasa(:,1:nocc),NTOvec) 1647if (nocc>10) then 1648 write(*,*) "The highest 10 eigenvalues of occupied NTOs:" 1649 write(*,"(5f12.6)") MOene(nocc-9:nocc) 1650else 1651 write(*,*) "Eigenvalues of occupied NTOs:" 1652 write(*,"(5f12.6)") MOene(1:nocc) 1653end if 1654deallocate(TT,NTOvec,NTOval) 1655!Virtual part 1656allocate(TT(nvir,nvir),NTOvec(nvir,nvir),NTOval(nvir)) 1657TT=matmul(transpose(T_MO),T_MO) 1658call diagsymat(TT,NTOvec,NTOval,ierror) 1659NTOval=NTOval*NTOvalcoeff 1660MOene(nocc+1:nbasis)=NTOval 1661CObasa(:,nocc+1:nbasis)=matmul(CObasa(:,nocc+1:nbasis),NTOvec) 1662do itmp=1,int(nvir/2D0) !Exchange array, so that the sequence will be high->low rather than the default low->high 1663 i=nocc+itmp 1664 j=nbasis+1-itmp 1665 tmpval=MOene(i) 1666 MOene(i)=MOene(j) 1667 MOene(j)=tmpval 1668 tmparr=CObasa(:,i) 1669 CObasa(:,i)=CObasa(:,j) 1670 CObasa(:,j)=tmparr 1671end do 1672write(*,*) 1673if (nvir>10) then 1674 write(*,*) "The highest 10 eigenvalues of virtual NTOs:" 1675 write(*,"(5f12.6)") MOene(nocc+1:nocc+10) 1676else 1677 write(*,*) "Eigenvalues of virtual NTOs:" 1678 write(*,"(5f12.6)") MOene(nocc+1:nbasis) 1679end if 1680deallocate(TT,NTOvec,NTOval,T_MO) 1681 1682!*** Beta part 1683if (allocated(CObasb)) then 1684 write(*,*) 1685 write(*,*) "Result of Beta part:" 1686 nocc=nint(nbelec) 1687 nvir=nbasis-nocc 1688 allocate(T_MO(nocc,nvir)) 1689 T_MO=0 1690 do iexcitorb=1,nexcitorb 1691 iocc=orbleft(iexcitorb) 1692 ivir=orbright(iexcitorb)-nocc 1693 if (iocc>nbasis) then !Only process transition of Beta orbitals 1694 iocc=iocc-nbasis 1695 ivir=ivir-nbasis 1696 if (excdir(iexcitorb)==1) T_MO(iocc,ivir)=T_MO(iocc,ivir)+exccoeff(iexcitorb) 1697 end if 1698 end do 1699 !Occupied part 1700 allocate(TT(nocc,nocc),NTOvec(nocc,nocc),NTOval(nocc)) 1701 TT=matmul(T_MO,transpose(T_MO)) 1702 call diagsymat(TT,NTOvec,NTOval,ierror) 1703 NTOval=NTOval*NTOvalcoeff 1704 MOene(nbasis+1:nbasis+nocc)=NTOval !By default, the diagsymat gives result from low to high 1705 CObasb(:,1:nocc)=matmul(CObasb(:,1:nocc),NTOvec) 1706 if (nocc>10) then 1707 write(*,*) "The highest 10 eigenvalues of occupied NTOs:" 1708 write(*,"(5f12.6)") MOene(nbasis+nocc-9:nbasis+nocc) 1709 else 1710 write(*,*) "Eigenvalues of occupied NTOs:" 1711 write(*,"(5f12.6)") MOene(nbasis+1:nbasis+nocc) 1712 end if 1713 deallocate(TT,NTOvec,NTOval) 1714 !Virtual part 1715 allocate(TT(nvir,nvir),NTOvec(nvir,nvir),NTOval(nvir)) 1716 TT=matmul(transpose(T_MO),T_MO) 1717 call diagsymat(TT,NTOvec,NTOval,ierror) 1718 NTOval=NTOval*NTOvalcoeff 1719 MOene(nbasis+nocc+1:nbasis+nbasis)=NTOval 1720 CObasb(:,nocc+1:nbasis)=matmul(CObasb(:,nocc+1:nbasis),NTOvec) 1721 do itmp=1,int(nvir/2D0) !Exchange array, so that the sequence will be high->low rather than the default low->high 1722 i=nocc+itmp 1723 j=nbasis+1-itmp 1724 tmpval=MOene(nbasis+i) 1725 MOene(nbasis+i)=MOene(nbasis+j) 1726 MOene(nbasis+j)=tmpval 1727 tmparr=CObasb(:,i) 1728 CObasb(:,i)=CObasb(:,j) 1729 CObasb(:,j)=tmparr 1730 end do 1731 write(*,*) 1732 if (nvir>10) then 1733 write(*,*) "The highest 10 eigenvalues of virtual NTOs:" 1734 write(*,"(5f12.6)") MOene(nbasis+nocc+1:nbasis+nocc+10) 1735 else 1736 write(*,*) "Eigenvalues of virtual NTOs:" 1737 write(*,"(5f12.6)") MOene(nbasis+nocc+1:nbasis+nbasis) 1738 end if 1739 deallocate(TT,NTOvec,NTOval,T_MO) 1740end if 1741write(*,*) 1742write(*,*) "0 Return" 1743write(*,*) "1 Output NTO orbitals to .molden file" 1744write(*,*) "2 Output NTO orbitals to .fch file" 1745read(*,*) iselNTO 1746if (iselNTO==1) then 1747 write(*,*) "Input the file path to output, e.g. C:\S1.molden" 1748 read(*,"(a)") c200tmp 1749 call outmolden(c200tmp,10) 1750 write(*,*) "Now you can load the newly generated .molden file to visualize NTOs" 1751else if (iselNTO==2) then 1752 write(*,*) "Input the file path to output, e.g. C:\S1.fch" 1753 read(*,"(a)") c200tmp 1754 call outfch(c200tmp,10) 1755 write(*,*) "Now you can load the newly generated .fch file to visualize NTOs" 1756end if 1757write(*,*) 1758write(*,*) "Reloading the initial file to recover status..." 1759call dealloall 1760call readinfile(firstfilename,1) 1761write(*,*) "Loading finished!" 1762write(*,*) 1763end subroutine 1764 1765 1766 1767 1768!!---------- Calculate all transition dipole moments between all excited states 1769!Note that this function cannot borrow the code in hetransdipdens for loading data, because this function & 1770!need transition information of all states as well as excitation energies 1771subroutine exctransdip 1772use defvar 1773use util 1774use function 1775implicit real*8 (a-h,o-z) 1776integer itype 1777!The last index in each arrays is the state index 1778integer,allocatable :: excdir(:,:) !excdir=1 means ->, =2 means <- 1779integer,allocatable :: orbleft(:,:),orbright(:,:) !denote the actual MO (have already considered alpha/beta problem) at the left/right side in the excitation data 1780real*8,allocatable :: exccoeff(:,:) !Coefficient of an orbital pair transition 1781real*8,allocatable :: excenergy(:) !Excitation energies 1782integer,allocatable :: excmulti(:) !Multiplicity of the states 1783integer,allocatable :: nexcitorb(:) !The number of MO pairs in the states 1784character c80tmp*80,transmodestr*200,leftstr*80,rightstr*80 1785character,save :: excitfilename*200=" " 1786real*8,allocatable :: GTFdipint(:,:) !Dipole moment integral between GTFs, use compressed index. The first index is x,y,z 1787real*8,allocatable :: MOdipint(:,:,:) !Dipole moment integral between all MOs. The first index is x,y,z 1788real*8,allocatable :: tdvecmat(:,:,:) 1789real*8 tdvec(3),grounddip(3),tmpvec(3) 1790 1791! excitfilename="c:\gtest\acetic_acid.out" 1792if (excitfilename==" ") then 1793 write(*,*) "Input the path of the Gaussian/ORCA output file or plain text file" 1794 write(*,*) "e.g. c:\lovelive\sunshine\yosoro.out" 1795 do while(.true.) 1796 read(*,"(a)") excitfilename 1797 inquire(file=excitfilename,exist=alive) 1798 if (alive) exit 1799 write(*,*) "Cannot find this file, input again" 1800 end do 1801else 1802 write(*,"(' Loading ',a)") trim(excitfilename) 1803end if 1804open(10,file=excitfilename,status="old") 1805call loclabel(10,"Gaussian",igauout,maxline=50) 1806call loclabel(10,"O R C A",iORCAout,maxline=50) 1807rewind(10) 1808 1809!Determine the number of excited states, so that proper size of arrays can be allocated 1810nstates=0 1811if (igauout==1) then !Gaussian output file 1812 write(*,*) "This is a Gaussian output file" 1813 call loclabel(10,"Excitation energies and oscillator strengths:") 1814 do while(.true.) 1815 call loclabel(10,"Excited State",ifound,0) 1816 if (ifound==1) then 1817 nstates=nstates+1 1818 read(10,*) 1819 else 1820 exit 1821 end if 1822 end do 1823else if (iORCAout==1) then !ORCA output file 1824 write(*,*) "This is an ORCA output file" 1825 call loclabel(10,"Number of roots to be determined",ifound) 1826 read(10,"(50x,i7)") nstates 1827else !Plain text file 1828 do while(.true.) 1829 call loclabel(10,"Excited State",ifound,0) 1830 if (ifound==1) then 1831 read(10,*) 1832 nstates=nstates+1 1833 else 1834 exit 1835 end if 1836 end do 1837end if 1838write(*,"(' There are',i5,' transitions, loading...')") nstates 1839 1840allocate(excenergy(nstates),excmulti(nstates),nexcitorb(nstates),tdvecmat(3,nstates,nstates)) 1841nexcitorb=0 1842 1843!Load excitation energy, multiplicity, the number of MO pairs of each excited state 1844if (igauout==1) then !Gaussian output file 1845 call loclabel(10,"Excitation energies and oscillator strengths:") 1846 do iexc=1,nstates 1847 call loclabel(10,"Excited State",ifound,0) 1848 read(10,"(a)") transmodestr 1849 excmulti(iexc)=0 !Multiplicity of the excited state 1850 if (index(transmodestr,"Singlet")/=0) excmulti(iexc)=1 1851 if (index(transmodestr,"Triplet")/=0) excmulti(iexc)=3 1852 do i=10,70 1853 if (transmodestr(i:i+1)=="eV") exit 1854 end do 1855 read(transmodestr(i-10:i-1),*) excenergy(iexc) 1856 !Count how many orbital pairs are involved in this transition mode 1857 do while(.true.) 1858 read(10,"(a)",iostat=ierror) c80tmp 1859 if (index(c80tmp,'<-')==0.and.index(c80tmp,'->')==0) exit 1860 nexcitorb(iexc)=nexcitorb(iexc)+1 1861 end do 1862 end do 1863else if (iORCAout==1) then !ORCA output file 1864 excmulti=1 !Multiplicity of all excited states are assumed to be singlet 1865 if (wfntype==0.or.wfntype==3) then 1866 call loclabel(10,"Generation of triplets") 1867 read(10,"(a)") c80tmp 1868 if (index(c80tmp," on ")/=0) then 1869 write(*,*) "Load which kind of excited states?" 1870 write(*,*) "1: Singlet 3: Triplet" 1871 read(*,*) iexcmulti 1872 excmulti=iexcmulti 1873 end if 1874 end if 1875 call loclabel(10,"the weight of the individual excitations are printed") 1876 if (iexcmulti==3) then !When triplets=on, ORCA calculate both singlet and triplet excited state, now move to the latter 1877 read(10,*) 1878 call loclabel(10,"the weight of the individual excitations are printed",ifound,0) 1879 end if 1880 do iexc=1,nstates 1881 call loclabel(10,"STATE",ifound,0) 1882 read(10,"(a)") transmodestr 1883 do i=10,70 1884 if (transmodestr(i:i+1)=="eV") exit 1885 end do 1886 read(transmodestr(i-10:i-1),*) excenergy(iexc) 1887 !Count how many orbital pairs are involved in this transition mode 1888 do while(.true.) 1889 read(10,"(a)",iostat=ierror) c80tmp 1890 if (index(c80tmp,'<-')==0.and.index(c80tmp,'->')==0) exit 1891 nexcitorb(iexc)=nexcitorb(iexc)+1 1892 end do 1893 end do 1894else 1895 rewind(10) 1896 do iexc=1,nstates 1897 call loclabel(10,"Excited State",ifound,0) 1898 read(10,*) c80tmp,c80tmp,inouse,excmulti(iexc),excenergy(iexc) 1899 !Count how many orbital pairs are involved in this transition mode 1900 do while(.true.) 1901 read(10,"(a)",iostat=ierror) c80tmp 1902 if ((index(c80tmp,'<-')==0.and.index(c80tmp,'->')==0).or.ierror/=0) exit 1903 nexcitorb(iexc)=nexcitorb(iexc)+1 1904 end do 1905 end do 1906end if 1907maxpair=maxval(nexcitorb) 1908 1909allocate(excdir(maxpair,nstates),orbleft(maxpair,nstates),orbright(maxpair,nstates),exccoeff(maxpair,nstates)) 1910 1911!Load MO transition coefficients and direction 1912if (igauout==1) then 1913 call loclabel(10,"Excitation energies and oscillator strengths:") 1914else if (iORCAout==1) then 1915 call loclabel(10,"the weight of the individual excitations are printed") 1916 if (iexcmulti==3) then !When triplets=on, ORCA calculate both singlet and triplet excited state, now move to the latter 1917 read(10,*) 1918 call loclabel(10,"the weight of the individual excitations are printed",ifound,0) 1919 end if 1920else 1921 rewind(10) 1922end if 1923!Notice that for unrestricted case, A and B are separately recorded in input file, & 1924!however after loading, they are combined as single index, namely if orbital index is larger than nbasis, then it is B, else A 1925if (iORCAout==1) then !ORCA output file 1926 !Worthnotingly, in at least ORCA 4.0, de-excitation is not separately outputted as <-, but combined into -> 1927 !Here we still determine <-, because hopefully Neese may change the convention of ORCA output in the future... 1928 do iexc=1,nstates 1929 call loclabel(10,"STATE",ifound,0) 1930 read(10,*) 1931 do itmp=1,nexcitorb(iexc) 1932 read(10,"(a)") c80tmp 1933 excdir(itmp,iexc)=-1 1934 if (index(c80tmp,'->')/=0) excdir(itmp,iexc)=1 !means -> 1935 if (index(c80tmp,'<-')/=0) excdir(itmp,iexc)=2 !means <- 1936! write(*,*) trim(c80tmp) 1937 do isign=1,80 !Find position of <- or -> 1938 if (c80tmp(isign:isign)=='-'.or.c80tmp(isign:isign)=='<') exit 1939 end do 1940 !Process left side of <- or -> 1941 read(c80tmp(:isign-1),"(a)") leftstr 1942 read(leftstr(:len_trim(leftstr)-1),*) orbleft(itmp,iexc) 1943 orbleft(itmp,iexc)=orbleft(itmp,iexc)+1 !ORCA counts orbital from 0 rather than 1!!! 1944 if (index(leftstr,'b')/=0) orbleft(itmp,iexc)=orbleft(itmp,iexc)+nbasis 1945 !Process right side of <- or -> 1946 read(c80tmp(isign+2:),*) rightstr 1947 read(rightstr(:len_trim(rightstr)-1),*) orbright(itmp,iexc) 1948 orbright(itmp,iexc)=orbright(itmp,iexc)+1 1949 if (index(rightstr,'b')/=0) orbright(itmp,iexc)=orbright(itmp,iexc)+nbasis 1950 iTDA=index(c80tmp,'c=') 1951 if (iTDA/=0) then !CIS, TDA task, configuration coefficients are presented 1952 read(c80tmp(iTDA+2:iTDA+13),*) exccoeff(itmp,iexc) 1953 else !TD task, configuration coefficients are not presented. Contribution of i->a and i<-a are summed up and outputted as i->a 1954 if (iexc==1.and.itmp==1) then 1955 write(*,"(a)") " Warning: For TD task, ORCA does not print configuration coefficients but only print corresponding contributions of each orbital pair, & 1956 in this case Multiwfn determines configuration coefficients simply as square root of contribution values. However, this treatment is & 1957 evidently inappropriate and the result is nonsense when de-excitation is significant (In this situation you have to use TDA-DFT instead)" 1958 write(*,*) "If you really want to proceed, press ENTER to continue" 1959 read(*,*) 1960 end if 1961 read(c80tmp(23:32),*) tmpval 1962 if (tmpval<0) excdir(itmp,iexc)=2 !Negative contribution is assumed to be de-excitation (of course this is not strict since -> and <- have been combined together) 1963 exccoeff(itmp,iexc)=dsqrt(abs(tmpval)) 1964 end if 1965 !Although for closed-shell ground state, ORCA still outputs coefficients as normalization to 100%, & 1966 !However, in order to follow the Gaussian convention, we change the coefficient as normalization to 50% 1967 if (wfntype==0.or.wfntype==3) exccoeff(itmp,iexc)=exccoeff(itmp,iexc)/dsqrt(2D0) 1968 end do 1969 end do 1970else !Gaussian output or plain text file 1971 do iexc=1,nstates 1972 call loclabel(10,"Excited State",ifound,0) 1973 read(10,*) 1974 do itmp=1,nexcitorb(iexc) 1975 read(10,"(a)") c80tmp 1976 excdir(itmp,iexc)=1 !means -> 1977 if (index(c80tmp,'<-')/=0) excdir(itmp,iexc)=2 !means <- 1978 do isign=1,80 !Find position of <- or -> 1979 if (c80tmp(isign:isign)=='-'.or.c80tmp(isign:isign)=='<') exit 1980 end do 1981 !Process left side of <- or -> 1982 read(c80tmp(:isign-1),"(a)") leftstr 1983 ilefttype=0 !close 1984 if (index(leftstr,'A')/=0) ilefttype=1 !Alpha 1985 if (index(leftstr,'B')/=0) ilefttype=2 !Beta 1986 if (ilefttype==0) then 1987 read(leftstr,*) orbleft(itmp,iexc) 1988 else 1989 read(leftstr(:len_trim(leftstr)-1),*) orbleft(itmp,iexc) 1990 if (ilefttype==2) orbleft(itmp,iexc)=orbleft(itmp,iexc)+nbasis 1991 end if 1992 !Process right side of <- or -> 1993 read(c80tmp(isign+2:),"(a)") rightstr 1994 irighttype=0 !close 1995 if (index(rightstr,'A')/=0) irighttype=1 !Alpha 1996 if (index(rightstr,'B')/=0) irighttype=2 !Beta 1997 if (irighttype==0) then 1998 read(rightstr,*) orbright(itmp,iexc),exccoeff(itmp,iexc) 1999 else 2000 do isplit=1,80 2001 if (rightstr(isplit:isplit)=='A'.or.rightstr(isplit:isplit)=='B') exit 2002 end do 2003 read(rightstr(:isplit-1),*) orbright(itmp,iexc) 2004 read(rightstr(isplit+1:),*) exccoeff(itmp,iexc) 2005 if (irighttype==2) orbright(itmp,iexc)=orbright(itmp,iexc)+nbasis 2006 end if 2007 end do 2008 end do 2009end if 2010close(10) 2011 2012!Test sum of square of the coefficients 2013write(*,*) "Summary:" 2014write(*,*) "Exc.state# Exc.energy(eV) Multi. N_pairs Sum coeff.^2" 2015do iexc=1,nstates 2016 sumsqrexc=0 2017 sumsqrdeexc=0 2018 do itmp=1,nexcitorb(iexc) 2019 if (excdir(itmp,iexc)==1) sumsqrexc=sumsqrexc+exccoeff(itmp,iexc)**2 2020 if (excdir(itmp,iexc)==2) sumsqrdeexc=sumsqrdeexc-exccoeff(itmp,iexc)**2 2021 end do 2022 sumsqrall=sumsqrexc+sumsqrdeexc 2023 write(*,"(i8,f18.5,4x,i8,i11,f16.6)") iexc,excenergy(iexc),excmulti(iexc),nexcitorb(iexc),sumsqrall 2024end do 2025write(*,*) 2026write(*,*) "Please select the destination of the output:" 2027write(*,*) "1 Output transition dipole moments to screen" 2028write(*,*) "2 Output transition dipole moments to transdipmom.txt in current folder" 2029write(*,*) "3 Generate input file of SOS module of Multiwfn as SOS.txt in current folder" 2030read(*,*) isel 2031write(*,*) 2032write(*,*) "Stage 1: Calculating dipole moment integrals between all GTFs..." 2033nsize=nprims*(nprims+1)/2 2034allocate(GTFdipint(3,nsize)) 2035call genGTFDmat(GTFdipint,nsize) 2036 2037call walltime(iwalltime1) 2038write(*,*) "Stage 2: Calculating dipole moment integrals between all MOs..." 2039allocate(MOdipint(3,nmo,nmo)) 2040!MOdipint will record dipole moment integrals between all MOs, including all occ+vir alpha and occ+vir beta ones 2041iprog=0 2042nthreads=getNThreads() 2043!$OMP PARALLEL DO SHARED(MOdipint,MOdipintb,iprog) PRIVATE(imo,jmo,iGTF,jGTF,ides,tmpvec) schedule(dynamic) NUM_THREADS(nthreads) 2044do imo=1,nmo 2045 do jmo=imo,nmo 2046 tmpvec=0 2047 do iGTF=1,nprims 2048 do jGTF=1,nprims 2049 if (iGTF>=jGTF) then 2050 ides=iGTF*(iGTF-1)/2+jGTF 2051 else 2052 ides=jGTF*(jGTF-1)/2+iGTF 2053 end if 2054 tmpvec=tmpvec+co(imo,iGTF)*co(jmo,jGTF)*GTFdipint(:,ides) 2055 end do 2056 end do 2057 MOdipint(:,imo,jmo)=tmpvec 2058 end do 2059 if (nprims>300) then 2060!$OMP CRITICAL 2061 iprog=iprog+1 2062 write(*,"(' Finished:',i6,' /',i6)") iprog,nmo 2063!$OMP END CRITICAL 2064 end if 2065end do 2066!$OMP END PARALLEL DO 2067!Fill lower triangle part 2068do imo=1,nmo 2069 do jmo=imo+1,nmo 2070 MOdipint(:,jmo,imo)=MOdipint(:,imo,jmo) 2071 end do 2072end do 2073call walltime(iwalltime2) 2074write(*,"(' (Stage 2 took up wall clock time',i10,'s)')") iwalltime2-iwalltime1 2075 2076if (isel==1) then 2077 iout=6 2078else if (isel==2) then 2079 iout=10 2080 open(iout,file="transdipmom.txt",status="replace") 2081else if (isel==3) then 2082 iout=10 2083 open(iout,file="SOS.txt",status="replace") 2084end if 2085fac=1 2086if (wfntype==0.or.wfntype==3) fac=2 2087grounddip=0 2088do imo=1,nmo 2089 grounddip=grounddip+MOocc(imo)*MOdipint(:,imo,imo) 2090end do 2091 2092!Output index, excitation energies and transition dipole moment between ground state and excited states, the format is in line with SOS module 2093if (isel==3) then 2094 write(iout,*) nstates 2095 do i=1,nstates 2096 write(iout,"(i6,f12.6)") i,excenergy(i) 2097 end do 2098 do iexc=0,nstates 2099 if (iexc==0) then 2100 write(iout,"(2i6,3(1PE15.6))") 0,iexc,grounddip 2101 else 2102 tdvec=0 2103 do ipair=1,nexcitorb(iexc) 2104 imo=orbleft(ipair,iexc) 2105 lmo=orbright(ipair,iexc) 2106 wei=exccoeff(ipair,iexc) 2107 tdvec=tdvec+wei*MOdipint(:,imo,lmo) 2108 end do 2109 write(iout,"(2i6,3(1PE15.6))") 0,iexc,tdvec*fac 2110 end if 2111 end do 2112end if 2113 2114write(*,*) "Stage 3: Calculating transition dipole moment between excited states..." 2115call walltime(iwalltime1) 2116iprog=0 2117nthreads=getNThreads() 2118!$OMP PARALLEL DO SHARED(tdvecmat,iprog) PRIVATE(iexc,jexc,tdvec,imo,lmo,jmo,kmo,wei) schedule(dynamic) NUM_THREADS(nthreads) 2119do iexc=1,nstates 2120 do jexc=iexc,nstates 2121 tdvec=0 2122 do ipair=1,nexcitorb(iexc) 2123 imo=orbleft(ipair,iexc) 2124 lmo=orbright(ipair,iexc) 2125 do jpair=1,nexcitorb(jexc) 2126 jmo=orbleft(jpair,jexc) 2127 kmo=orbright(jpair,jexc) 2128 wei=exccoeff(ipair,iexc)*exccoeff(jpair,jexc) 2129 if (excdir(ipair,iexc)==excdir(jpair,jexc)) then !If don't apply this condition, for TD may be evidently deviate to actual result 2130 if (excdir(ipair,iexc)==1) then ! -> case 2131 if (imo==jmo.and.lmo/=kmo) then 2132 tdvec=tdvec+wei*MOdipint(:,lmo,kmo) 2133 else if (imo/=jmo.and.lmo==kmo) then 2134 tdvec=tdvec-wei*MOdipint(:,imo,jmo) 2135 else if (imo==jmo.and.lmo==kmo) then 2136 tdvec=tdvec+wei*(grounddip-MOdipint(:,imo,imo)+MOdipint(:,lmo,lmo)) 2137 end if 2138! else ! <- case. Frankly speaking, I don't know how to exactly deal with this circumstance.& 2139 !So I simply ignore them and I found this is the best way currently. The error must be very small, since de-excitation coefficients are often quite small 2140! if (imo==jmo.and.lmo/=kmo) then 2141! tdvec=tdvec-wei*MOdipint(:,lmo,kmo) 2142! else if (imo/=jmo.and.lmo==kmo) then 2143! tdvec=tdvec+wei*MOdipint(:,imo,jmo) 2144! else if (imo==jmo.and.lmo==kmo) then 2145! tdvec=tdvec+wei*(grounddip+MOdipint(:,imo,imo)-MOdipint(:,lmo,lmo)) 2146! end if 2147 end if 2148 end if 2149 end do 2150 end do 2151 tdvecmat(:,iexc,jexc)=tdvec*fac 2152 end do 2153 2154 if (nprims>300) then 2155!$OMP CRITICAL 2156 iprog=iprog+1 2157 write(*,"(' Finished:',i6,' /',i6)") iprog,nstates 2158!$OMP END CRITICAL 2159 end if 2160end do 2161!$OMP END PARALLEL DO 2162call walltime(iwalltime2) 2163write(*,"(' (Stage 3 took up wall clock time',i10,'s)',/)") iwalltime2-iwalltime1 2164 2165if (isel==1.or.isel==2) then 2166 write(iout,"(' Ground state dipole moment in X,Y,Z:',3f12.6,' a,u,',/)") grounddip 2167 write(iout,"(' Transition dipole moment between excited states (a.u.):')") 2168 write(iout,*) " i j X Y Z Diff.(eV) Oscil.str" 2169end if 2170 2171do iexc=1,nstates 2172 do jexc=iexc,nstates 2173 if (isel==1.or.isel==2) then 2174 enediff=abs(excenergy(jexc)-excenergy(iexc)) 2175 oscillstr=2D0/3D0*enediff/au2eV*sum(tdvecmat(:,iexc,jexc)**2) 2176 write(iout,"(2i6,3f14.7,2f12.5)") iexc,jexc,tdvecmat(:,iexc,jexc),enediff,oscillstr 2177 else if (isel==3) then 2178 write(iout,"(2i6,3(1PE15.6))") iexc,jexc,tdvecmat(:,iexc,jexc) 2179 end if 2180 end do 2181end do 2182 2183if (isel==2) then 2184 close(iout) 2185 write(*,*) "Done! The result has been outputted to transdipmom.txt in current folder" 2186else if (isel==3) then 2187 close(iout) 2188 write(*,*) "Done! The result has been outputted to SOS.txt in current folder" 2189end if 2190write(*,*) 2191end subroutine 2192 2193 2194 2195 2196!!!------------- Analyze charge transfer 2197!See J. Chem. Theory Comput., 7, 2498 2198subroutine CTanalyze 2199use defvar 2200implicit real*8(a-h,o-z) 2201real*8,allocatable :: Cpos(:,:,:),Cneg(:,:,:),tmpmat(:,:,:) 2202sumpos=0D0 2203sumneg=0D0 2204Rxpos=0D0 2205Rypos=0D0 2206Rzpos=0D0 2207Rxneg=0D0 2208Ryneg=0D0 2209Rzneg=0D0 2210do i=1,nx 2211 do j=1,ny 2212 do k=1,nz 2213 if (cubmat(i,j,k)>0) then 2214 sumpos=sumpos+cubmat(i,j,k) 2215 Rxpos=Rxpos+cubmat(i,j,k)*(orgx+(i-1)*dx) 2216 Rypos=Rypos+cubmat(i,j,k)*(orgy+(j-1)*dy) 2217 Rzpos=Rzpos+cubmat(i,j,k)*(orgz+(k-1)*dz) 2218 else 2219 sumneg=sumneg+cubmat(i,j,k) 2220 Rxneg=Rxneg+cubmat(i,j,k)*(orgx+(i-1)*dx) 2221 Ryneg=Ryneg+cubmat(i,j,k)*(orgy+(j-1)*dy) 2222 Rzneg=Rzneg+cubmat(i,j,k)*(orgz+(k-1)*dz) 2223 end if 2224 end do 2225 end do 2226end do 2227dvol=dx*dy*dz 2228sumpos=sumpos*dvol 2229sumneg=sumneg*dvol 2230Rxpos=Rxpos*dvol/sumpos 2231Rypos=Rypos*dvol/sumpos 2232Rzpos=Rzpos*dvol/sumpos 2233Rxneg=Rxneg*dvol/sumneg 2234Ryneg=Ryneg*dvol/sumneg 2235Rzneg=Rzneg*dvol/sumneg 2236disx=abs(Rxpos-Rxneg) 2237disy=abs(Rypos-Ryneg) 2238disz=abs(Rzpos-Rzneg) 2239disnorm=dsqrt(disx**2+disy**2+disz**2) 2240write(*,"(' Transferred charge (positive and negative parts):',2f8.3)") sumpos,sumneg 2241write(*,"(' Barycenter of positive part in x,y,z (Angstrom):',3f8.3)") Rxpos*b2a,Rypos*b2a,Rzpos*b2a 2242write(*,"(' Barycenter of negative part in x,y,z (Angstrom):',3f8.3)") Rxneg*b2a,Ryneg*b2a,Rzneg*b2a 2243write(*,"(' Distance of CT in x,y,z (Angstrom):',3f8.3,' Norm:',f8.3)") disx*b2a,disy*b2a,disz*b2a,disnorm*b2a 2244dipx=-(Rxpos-Rxneg)*sumpos 2245dipy=-(Rypos-Ryneg)*sumpos 2246dipz=-(Rzpos-Rzneg)*sumpos 2247dipnorm=disnorm*sumpos 2248write(*,"(' Dipole moment variation (a.u.) :',3f8.3,' Norm:',f8.3)") dipx,dipy,dipz,dipnorm 2249write(*,"(' Dipole moment variation (Debye):',3f8.3,' Norm:',f8.3)") dipx*au2debye,dipy*au2debye,dipz*au2debye,dipnorm*au2debye 2250 2251sigxpos=0D0 2252sigypos=0D0 2253sigzpos=0D0 2254sigxneg=0D0 2255sigyneg=0D0 2256sigzneg=0D0 2257do i=1,nx 2258 do j=1,ny 2259 do k=1,nz 2260 if (cubmat(i,j,k)>0) then 2261 sigxpos=sigxpos+cubmat(i,j,k)*((orgx+(i-1)*dx)-Rxpos)**2 2262 sigypos=sigypos+cubmat(i,j,k)*((orgy+(j-1)*dy)-Rypos)**2 2263 sigzpos=sigzpos+cubmat(i,j,k)*((orgz+(k-1)*dz)-Rzpos)**2 2264 else 2265 sigxneg=sigxneg+cubmat(i,j,k)*((orgx+(i-1)*dx)-Rxneg)**2 2266 sigyneg=sigyneg+cubmat(i,j,k)*((orgy+(j-1)*dy)-Ryneg)**2 2267 sigzneg=sigzneg+cubmat(i,j,k)*((orgz+(k-1)*dz)-Rzneg)**2 2268 end if 2269 end do 2270 end do 2271end do 2272sigxpos=dsqrt(sigxpos/(sumpos/dvol)) 2273sigypos=dsqrt(sigypos/(sumpos/dvol)) 2274sigzpos=dsqrt(sigzpos/(sumpos/dvol)) 2275signormpos=dsqrt(sigxpos**2+sigypos**2+sigzpos**2) 2276sigxneg=dsqrt(sigxneg/(sumneg/dvol)) 2277sigyneg=dsqrt(sigyneg/(sumneg/dvol)) 2278sigzneg=dsqrt(sigzneg/(sumneg/dvol)) 2279signormneg=dsqrt(sigxneg**2+sigyneg**2+sigzneg**2) 2280write(*,"(' RMSD of positive part in x,y,z (Angstrom):',3f8.3,' Tot:',f7.3)") sigxpos*b2a,sigypos*b2a,sigzpos*b2a,signormpos*b2a 2281write(*,"(' RMSD of negative part in x,y,z (Angstrom):',3f8.3,' Tot:',f7.3)") sigxneg*b2a,sigyneg*b2a,sigzneg*b2a,signormneg*b2a 2282delta_sigCT=dsqrt((sigxpos-sigxneg)**2+(sigypos-sigyneg)**2+(sigzpos-sigzneg)**2) 2283write(*,"(' Difference between RMSD of positive and negative parts (Angstrom):',f8.3)") delta_sigCT*b2a 2284Hx=(sigxpos+sigxneg)/2D0 2285Hy=(sigypos+sigyneg)/2D0 2286Hz=(sigzpos+sigzneg)/2D0 2287Hnorm=dsqrt(Hx**2+Hy**2+Hz**2) 2288write(*,"(' H index in x,y,z (Angstrom):',3f8.3,' Norm:',f8.3)") Hx*b2a,Hy*b2a,Hz*b2a,Hnorm*b2a 2289tx=disx-Hx 2290ty=disy-Hy 2291tz=disz-Hz 2292tnorm=dsqrt(tx**2+ty**2+tz**2) 2293write(*,"(' t index in x,y,z (Angstrom):',3f8.3,' Norm:',f8.3)") tx*b2a,ty*b2a,tz*b2a,tnorm*b2a 2294 2295allocate(Cpos(nx,ny,nz),Cneg(nx,ny,nz),tmpmat(nx,ny,nz)) 2296do i=1,nx 2297 do j=1,ny 2298 do k=1,nz 2299 rnowx=orgx+(i-1)*dx 2300 rnowy=orgy+(j-1)*dy 2301 rnowz=orgz+(k-1)*dz 2302 Cpos(i,j,k)=exp( -(rnowx-Rxpos)**2/(2*sigxpos**2) -(rnowy-Rypos)**2/(2*sigypos**2) -(rnowz-Rzpos)**2/(2*sigzpos**2)) 2303 Cneg(i,j,k)=exp( -(rnowx-Rxneg)**2/(2*sigxneg**2) -(rnowy-Ryneg)**2/(2*sigyneg**2) -(rnowz-Rzneg)**2/(2*sigzneg**2)) 2304 end do 2305 end do 2306end do 2307Cpos=Cpos*sumpos/(sum(Cpos)*dvol) 2308Cneg=Cneg*sumneg/(sum(Cneg)*dvol) 2309 2310write(*,"(' Overlap integral between C+ and C-:',f12.6)") sum(dsqrt(Cpos/sumpos)*dsqrt(Cneg/sumneg))*dvol 2311 2312sur_value=0.001 2313do while(.true.) 2314 write(*,*) 2315 write(*,*) "0 Return" 2316 write(*,*) "1 Show isosurface of C+ and C- functions simultaneously" 2317 write(*,*) "2 Export C+ and C- functions to cube file in current folder" 2318 read(*,*) isel 2319 2320 if (isel==0) then 2321 return 2322 else if (isel==1) then 2323 isosursec=1 2324 clrRcub1sameold=clrRcub1same !Backup previous color setting 2325 clrGcub1sameold=clrGcub1same 2326 clrBcub1sameold=clrBcub1same 2327 clrRcub2oppoold=clrRcub2oppo 2328 clrGcub2oppoold=clrGcub2oppo 2329 clrBcub2oppoold=clrBcub2oppo 2330 clrRcub2samemeshptold=clrRcub2samemeshpt 2331 clrGcub2samemeshptold=clrGcub2samemeshpt 2332 clrBcub2samemeshptold=clrBcub2samemeshpt 2333 clrRcub2oppomeshptold=clrRcub2oppomeshpt 2334 clrGcub2oppomeshptold=clrGcub2oppomeshpt 2335 clrBcub2oppomeshptold=clrBcub2oppomeshpt 2336 clrRcub1same=0.3D0 !Set color to that Cpos is green, Cneg is blue. (Cpos/Cneg function is positive/negative everywhere due to sumpos and sumneg) 2337 clrGcub1same=0.75D0 2338 clrBcub1same=0.3D0 2339 clrRcub2oppo=0.3D0 2340 clrGcub2oppo=0.45D0 2341 clrBcub2oppo=0.9D0 2342 clrRcub2samemeshpt=0.3D0 2343 clrGcub2samemeshpt=0.75D0 2344 clrBcub2samemeshpt=0.3D0 2345 clrRcub2oppomeshpt=0.3D0 2346 clrGcub2oppomeshpt=0.45D0 2347 clrBcub2oppomeshpt=0.9D0 2348 if (allocated(cubmattmp)) deallocate(cubmattmp) 2349 allocate(cubmattmp(nx,ny,nz)) 2350 tmpmat=cubmat 2351 cubmat=Cpos 2352 cubmattmp=Cneg 2353 !Since drawisosurgui is mainly designed for plotting one isosurface, while here we use it to draw both cubmat and cubmattmp, so disable change of style to avoid troubles 2354 cubmat=tmpmat !Recover original data, namely density difference between two states 2355 deallocate(cubmattmp) 2356 clrRcub1same=clrRcub1sameold !Recover previous color setting 2357 clrGcub1same=clrGcub1sameold 2358 clrBcub1same=clrBcub1sameold 2359 clrRcub2oppo=clrRcub2oppoold 2360 clrGcub2oppo=clrGcub2oppoold 2361 clrBcub2oppo=clrBcub2oppoold 2362 clrRcub2samemeshpt=clrRcub2samemeshptold 2363 clrGcub2samemeshpt=clrGcub2samemeshptold 2364 clrBcub2samemeshpt=clrBcub2samemeshptold 2365 clrRcub2oppomeshpt=clrRcub2oppomeshptold 2366 clrGcub2oppomeshpt=clrGcub2oppomeshptold 2367 clrBcub2oppomeshpt=clrBcub2oppomeshptold 2368 else if (isel==2) then 2369 open(10,file="Cpos.cub",status="replace") 2370 call outcube(Cpos,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 2371 close(10) 2372 write(*,"('C+ function has been outputted to ""Cpos.cub"" in current folder')") 2373 open(10,file="Cneg.cub",status="replace") 2374 call outcube(Cneg,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 2375 close(10) 2376 write(*,"('C- function has been outputted to ""Cneg.cub"" in current folder')") 2377 end if 2378end do 2379end subroutine 2380 2381 2382 2383 2384 2385!-------------------------------------------------------------------------- 2386!----------- Plot transition density matrix as color-filled map ----------- 2387!-------------------------------------------------------------------------- 2388subroutine plottransdensmat 2389use defvar 2390use util 2391implicit real*8 (a-h,o-z) 2392character tdmatfilename*200 2393real*8 tdmatbas(nbasis,nbasis),tdmatatmtmp(ncenter,ncenter),tdmatatm(ncenter,ncenter) 2394real*8,allocatable :: tdmatatmnoh(:,:) 2395!Load density transition matrix in basis expansion 2396tdmatbas=0D0 2397write(*,"(a)") " Input the path of the Gaussian output file or plain text file containing transition density matrix, e.g. c:\a.out" 2398do while(.true.) 2399 read(*,"(a)") tdmatfilename 2400 inquire(file=tdmatfilename,exist=alive) 2401 if (alive) exit 2402 write(*,*) "Error: Cannot find the file, please input again" 2403end do 2404 2405open(10,file=tdmatfilename,status="old") 2406if (index(tdmatfilename,"AAtrdip.txt")/=0) then !Directly load and use atom-atom matrix (commonly generated by option 5 of hole-electron module) 2407 write(*,*) "Loading atom-atom contribution matrix..." 2408 call readmatgau(10,tdmatatm,0,"f14.8",6,5,1) 2409else !Load basis-basis matrix and then condense to atom-atom matrix 2410 call loclabel(10,"Gaussian",igauout,maxline=100) 2411 rewind(10) 2412 if (igauout==1) then !Gaussian output file 2413 write(*,*) "This is a Gaussian output file" 2414 call loclabel(10,"Alpha Density Matrix:",ifound) 2415 if (ifound==1) then !Open-shell 2416 write(*,*) "Use which type of transition density matrix?" 2417 write(*,*) "1=Alpha 2=Beta" 2418 read(*,*) iTDMtype 2419 write(*,*) "Loading transition density matrix..." 2420 if (iTDMtype==1) then 2421 call readmatgau(10,tdmatbas,1,"f10.5",21,5,1) 2422 else if (iTDMtype==2) then 2423 call loclabel(10,"Beta Density Matrix:",ifound) 2424 call readmatgau(10,tdmatbas,1,"f10.5",21,5,1) 2425 end if 2426 else !Close-shell 2427 call loclabel(10,"Density Matrix:",ifound) 2428 if (ifound==0) call loclabel(10,"DENSITY MATRIX.",ifound) 2429 if (ifound==0) then 2430 write(*,"(a,/)") "Error: Cannot found transition density matrix information from the Gaussian output file, please check if iop(6/8=3) has been specified" 2431 return 2432 end if 2433 write(*,*) "Loading transition density matrix..." 2434 call readmatgau(10,tdmatbas,1,"f10.5",21,5,1) 2435 end if 2436 else !Plain text file with transition density matrix in basis function 2437 write(*,*) "Loading transition density matrix..." 2438 call readmatgau(10,tdmatbas,0,"f14.8",6,5,1) 2439 end if 2440 !Use tdmatbas to construct tdmatatm 2441 tdmatatm=0D0 2442 do iatm=1,ncenter 2443 do jatm=1,ncenter 2444 do ibas=basstart(iatm),basend(iatm) 2445 do jbas=basstart(jatm),basend(jatm) 2446 tdmatatm(iatm,jatm)=tdmatatm(iatm,jatm)+tdmatbas(ibas,jbas)**2 2447 ! tdmatatm(iatm,jatm)=tdmatatm(iatm,jatm)+tdmatbas(ibas,jbas) 2448 ! tdmatatm(iatm,jatm)=tdmatatm(iatm,jatm)+abs(tdmatbas(ibas,jbas)) 2449 end do 2450 end do 2451 end do 2452 end do 2453end if 2454close(10) 2455 2456tdmatatmtmp=tdmatatm 2457!Contract tdmatatm to tdmatatmnoh, tdmatatmtmp is used as a intermediate 2458!tdmatatmnoh is the tdmatatm without hydrogens 2459ncenreal=count(a(:)%index/=1) 2460write(*,"(' The number of non-hydrogen atoms:',i10)") ncenreal 2461allocate(tdmatatmnoh(ncenreal,ncenreal)) 2462itmp=0 2463do iatm=1,ncenter 2464 if (a(iatm)%index/=1) then 2465 itmp=itmp+1 2466 tdmatatmtmp(:,itmp)=tdmatatmtmp(:,iatm) 2467 tdmatatmtmp(itmp,:)=tdmatatmtmp(iatm,:) 2468 end if 2469end do 2470tdmatatmnoh(:,:)=tdmatatmtmp(1:ncenreal,1:ncenreal) 2471 2472ifhydrogen=0 2473clrlimlow=minval(tdmatatmnoh) 2474clrlimhigh=maxval(tdmatatmnoh) 2475ninterpo=10 2476nstepsize=5 2477ifnormsum=0 2478! clrlimlow=minval(tdmatbas) 2479! clrlimhigh=maxval(tdmatbas) 2480! write(*,*) clrlimlow,clrlimhigh 2481facnorm=1D0 !Normalization factor, default is 1, namely don't do normalization 2482 2483do while(.true.) 2484 write(*,*) 2485 write(*,*) "0 Return" 2486 write(*,*) "1 Show transition density matrix map" 2487 write(*,*) "2 Save transition density matrix map to a graphical file in current folder" 2488 write(*,*) "3 Export data to plain text file" 2489 if (ifhydrogen==0) write(*,*) "4 Switch if take hydrogens into account, current: No" 2490 if (ifhydrogen==1) write(*,*) "4 Switch if take hydrogens into account, current: Yes" 2491 write(*,"(a,f7.4,a,f7.4)") " 5 Change lower and upper limit of color scale, current:",clrlimlow," to",clrlimhigh 2492 write(*,"(a,i3)") " 6 Set the number of interpolation steps between grid, current:",ninterpo 2493 write(*,"(a,i3)") " 7 Set stepsize between labels, current:",nstepsize 2494 if (ifnormsum==0) write(*,*) "8 Switch if normalized the sum of all elements to unity, current: No" 2495 if (ifnormsum==1) write(*,*) "8 Switch if normalized the sum of all elements to unity, current: Yes" 2496 read(*,*) isel 2497 2498 if (isel==0) then 2499 return 2500 else if (isel==1.or.isel==2) then 2501 if (isel==2) isavepic=1 2502 if (ifhydrogen==1) then 2503 nspc=(ncenter-1)/nlabstep 2504 else if (ifhydrogen==0) then 2505 nspc=(ncenreal-1)/nlabstep 2506 end if 2507 if (isel==2) then 2508 isavepic=0 2509 write(*,*) "Done, the picture has been saved to current folder with ""DISLIN"" prefix" 2510 end if 2511 else if (isel==3) then 2512 open(10,file="tdmat.txt",status="replace") 2513 if (ifhydrogen==0) then 2514 do iatm=1,ncenreal 2515 do jatm=1,ncenreal 2516 write(10,"(2i8,f12.6)") iatm,jatm,tdmatatmnoh(iatm,jatm)/facnorm 2517 end do 2518 end do 2519 else if (ifhydrogen==1) then 2520 do iatm=1,ncenter 2521 do jatm=1,ncenter 2522 write(10,"(2i8,f12.6)") iatm,jatm,tdmatatm(iatm,jatm)/facnorm 2523 end do 2524 end do 2525 end if 2526 close(10) 2527 write(*,*) "Done, the data have been exported to tdmat.txt in current folder" 2528 else if (isel==4) then 2529 if (ifhydrogen==0) then 2530 ifhydrogen=1 2531 clrlimlow=minval(tdmatatm) 2532 clrlimhigh=maxval(tdmatatm) 2533 else if (ifhydrogen==1) then 2534 ifhydrogen=0 2535 clrlimlow=minval(tdmatatmnoh) 2536 clrlimhigh=maxval(tdmatatmnoh) 2537 end if 2538 else if (isel==5) then 2539 write(*,*) "Input lower and upper limits, e.g. 0,1.5" 2540 read(*,*) clrlimlow,clrlimhigh 2541 else if (isel==6) then 2542 write(*,*) "Please input a number" 2543 write(*,"(a)") " Note: Larger value gives rise to more smooth graph, 1 means don't do interpolation" 2544 read(*,*) ninterpo 2545 else if (isel==7) then 2546 write(*,*) "Please input a number, e.g. 5" 2547 read(*,*) nstepsize 2548 else if (isel==8) then 2549 if (ifnormsum==0) then 2550 ifnormsum=1 2551 facnorm=sum(tdmatatm(:,:)) !The sum of all elements 2552 write(*,*) "The normalization factor is", facnorm 2553 else if (ifnormsum==1) then 2554 ifnormsum=0 2555 facnorm=1D0 2556 end if 2557 end if 2558end do 2559end subroutine 2560 2561 2562 2563 2564 2565!--------------------------------------------------------------------------------------- 2566!---------- Calculate interfragment charger transfer in electronic excitation ---------- 2567!--------------------------------------------------------------------------------------- 2568subroutine excfragCT(nexcitorb,orbleft,orbright,excdir,exccoeff) 2569use defvar 2570use util 2571implicit real*8 (a-h,o-z) 2572integer nexcitorb,orbleft(nexcitorb),orbright(nexcitorb),excdir(nexcitorb) 2573real*8 exccoeff(nexcitorb),CTmat(ncenter,ncenter),CTmattmp(ncenter,ncenter),atmcomp(ncenter,nmo) 2574character c2000tmp*2000 2575 2576inquire(file="orbcomp.txt",exist=alive) 2577if (alive) then !Load atomic contribution from orbcomp.txt, which may be outputted by option -4 of Hirshfeld/Becke composition analysis 2578 write(*,"(a)") " orbcomp.txt was found in current folder, now load atomic contribution to all orbitals from this file..." 2579 open(10,file="orbcomp.txt",status="old") 2580 do imo=1,nmo 2581 read(10,*) 2582 do iatm=1,ncenter 2583 read(10,*) inouse,atmcomp(iatm,imo) 2584 end do 2585 end do 2586 close(10) 2587 atmcomp=atmcomp/100 2588else 2589 write(*,*) "Calculating atomic contribution to all orbitals by SCPA method..." 2590 do imo=1,nmo 2591 if (MOtype(imo)==0.or.MOtype(imo)==1) then !Close-shell or alpha part of open-shell 2592 sumsqr=sum(cobasa(:,imo)**2) 2593 do iatm=1,ncenter 2594 atmcomp(iatm,imo)=sum(cobasa(basstart(iatm):basend(iatm),imo)**2)/sumsqr 2595 end do 2596 else !Beta part of open-shell 2597 iimo=imo-nbasis 2598 sumsqr=sum(cobasb(:,iimo)**2) 2599 do iatm=1,ncenter 2600 atmcomp(iatm,imo)=sum(cobasb(basstart(iatm):basend(iatm),iimo)**2)/sumsqr 2601 end do 2602 end if 2603 end do 2604end if 2605 2606write(*,*) "Constructing inter-atomic CT matrix..." 2607CTmat=0 2608fac=1 2609if (wfntype==0.or.wfntype==3) fac=2 !Since for closed-shell case, program only print one-half part 2610do iexc=1,nexcitorb 2611 imo=orbleft(iexc) 2612 jmo=orbright(iexc) 2613 !Calculate transferred electron from iatm(row) to jatm(column) 2614 do iatm=1,ncenter 2615 do jatm=1,ncenter 2616 CTmattmp(iatm,jatm)=atmcomp(iatm,imo)*atmcomp(jatm,jmo) 2617 end do 2618 end do 2619 if (excdir(iexc)/=1) CTmattmp=-CTmattmp 2620 CTmat=CTmat+CTmattmp*fac*exccoeff(iexc)**2 2621end do 2622write(*,*) "Preparation is finished!" 2623 2624do while(.true.) 2625 write(*,*) 2626 write(*,*) "Input atom list for fragment 1, e.g. 3,5-8,15-20" 2627 write(*,*) "Note: Input 0 can exit" 2628 read(*,"(a)") c2000tmp 2629 if (c2000tmp(1:1)=="0") return 2630 call str2arr(c2000tmp,nterm1) 2631 if (allocated(frag1)) deallocate(frag1) 2632 allocate(frag1(nterm1)) 2633 call str2arr(c2000tmp,nterm1,frag1) 2634 write(*,*) "Input atom list for fragment 2, e.g. 1,2,4,9-14" 2635 read(*,"(a)") c2000tmp 2636 call str2arr(c2000tmp,nterm2) 2637 if (allocated(frag2)) deallocate(frag2) 2638 allocate(frag2(nterm2)) 2639 call str2arr(c2000tmp,nterm2,frag2) 2640 2641 CTval1=0 2642 CTval2=0 2643 varpop1=0 2644 varpop2=0 2645 do idx=1,nterm1 2646 iatm=frag1(idx) 2647 do jdx=1,nterm2 2648 jatm=frag2(jdx) 2649 CTval1=CTval1+CTmat(iatm,jatm) 2650 CTval2=CTval2+CTmat(jatm,iatm) 2651 end do 2652 varpop1=varpop1+sum(CTmat(:,iatm))-sum(CTmat(iatm,:)) 2653 end do 2654 do jdx=1,nterm2 2655 jatm=frag2(jdx) 2656 varpop2=varpop2+sum(CTmat(:,jatm))-sum(CTmat(jatm,:)) 2657 end do 2658 write(*,"(a,f10.5)") " Electron transferred from fragment 1 to 2:",CTval1 2659 write(*,"(a,f10.5)") " Electron transferred from fragment 2 to 1:",CTval2 2660 write(*,"(a,f10.5)") " Electron net transferred from fragment 1 to 2:",CTval1-CTval2 2661 write(*,"(a,f10.5)") " Variation of electron population of fragment 1:",varpop1 2662 write(*,"(a,f10.5)") " Variation of electron population of fragment 2:",varpop2 2663end do 2664 2665end subroutine 2666