diff --git a/linear_response.f b/linear_response.f index b88df22..ce120f0 100644 --- a/linear_response.f +++ b/linear_response.f @@ -19,24 +19,24 @@ ! written by Marc de Wegifosse 2017-2019 - SUBROUTINE lresp(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) + SUBROUTINE lresp(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk integer ::maxconf,moci,no,nv integer ::iconf(maxconf,2) integer*8 ::lin8 - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) - - + + real*4 ::apb(nci*(nci+1)/2) real*4 ::amb(nci*(nci+1)/2) real*4, allocatable ::inv_amb(:) @@ -51,16 +51,16 @@ SUBROUTINE lresp(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) integer ::info integer, allocatable ::ipiv(:) real*4, allocatable ::work (:) - + integer ::ix,iy,iz - + real*4 ::start_time,end_time - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) mu(:,3)=zl(:) - + open(unit=101,file='wavelength',form='formatted',status='old') freq(1)=0.0 write(*,*) @@ -79,15 +79,15 @@ SUBROUTINE lresp(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) freq(i+1)=45.56335/freq(i+1) enddo close(101) - + allocate(inv_amb(nci*(nci+1)/2)) inv_amb=amb uplo='U' allocate(ipiv(1:nci),work(1:nci)) call ssptrf(uplo,nci,inv_amb,ipiv,info) - call ssptri(uplo,nci,inv_amb,ipiv,work,info) + call ssptri(uplo,nci,inv_amb,ipiv,work,info) deallocate(ipiv,work) - + allocate( XpY(nci,3)) Do ii=1, num_freq+1 @@ -96,16 +96,16 @@ SUBROUTINE lresp(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) call cpu_time(start_time) allocate(inv_resp(nci*(nci+1)/2)) inv_resp=apb-omega**2.0*inv_amb - + uplo='U' allocate(ipiv(1:nci),work(1:nci)) call ssptrf(uplo,nci,inv_resp,ipiv,info) - call ssptri(uplo,nci,inv_resp,ipiv,work,info) + call ssptri(uplo,nci,inv_resp,ipiv,work,info) deallocate(ipiv,work) - - - + + + !$omp parallel private(i,j,jk,io,iv,idum1,idum2,ij) !$omp& reduction (+:XpY) !$omp do @@ -139,9 +139,9 @@ SUBROUTINE lresp(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) io=iconf(j,1) iv=iconf(j,2) idum1=max(io,iv) - idum2=min(io,iv) + idum2=min(io,iv) ij=idum2+idum1*(idum1-1)/2 - alpha_xx=alpha_xx-(xl(ij)*XpY(j,1))*2.0 + alpha_xx=alpha_xx-(xl(ij)*XpY(j,1))*2.0 alpha_xy=alpha_xy-(xl(ij)*XpY(j,2))*2.0 alpha_xz=alpha_xz-(xl(ij)*XpY(j,3))*2.0 alpha_yy=alpha_yy-(yl(ij)*XpY(j,2))*2.0 @@ -165,19 +165,19 @@ SUBROUTINE lresp(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) enddo deallocate(XpY) deallocate(inv_amb) -111 format(A15,F20.6) +111 format(A15,F20.6) 1111 format(A22,2A20) -2222 format(A2,3F20.6) +2222 format(A2,3F20.6) 3333 format(A16,F7.1,A1,F7.1,A1) end subroutine lresp - - - - SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) + + + + SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk integer ::maxconf,moci,no,nv @@ -187,14 +187,14 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) integer, allocatable :: B_list(:,:) integer ::counter_B integer*8 ::lin8 - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) - - + + real*4 ::apb(nci*(nci+1)/2) real*4 ::amb(nci*(nci+1)/2) real*4, allocatable ::inv_amb(:) @@ -212,14 +212,14 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) integer ::info integer, allocatable ::ipiv(:) real*4, allocatable ::work (:) - + integer ::ix,iy,iz real*8 ::beta(3,3,3),A,B real*8 ::beta2_ZZZ, beta2_XZZ, beta_HRS, DR real*8 ::betaVEC_X, betaVEC_Y, betaVEC_Z - + real*4 ::start_time,end_time - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) @@ -245,15 +245,15 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) freq(i+1)=45.56335/freq(i+1) enddo close(101) - + allocate(inv_amb(nci*(nci+1)/2)) inv_amb=amb uplo='U' allocate(ipiv(1:nci),work(1:nci)) call ssptrf(uplo,nci,inv_amb,ipiv,info) - call ssptri(uplo,nci,inv_amb,ipiv,work,info) + call ssptri(uplo,nci,inv_amb,ipiv,work,info) deallocate(ipiv,work) - + allocate( XpY(nci,num_freq+1,2,3)) allocate( XmY(nci,num_freq+1,2,3)) allocate( X(nci,num_freq+1,2,3), @@ -266,20 +266,20 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) Do ii=1, num_freq+1 Do jj=1, 2 !2*omega omega=freq(ii)*jj - if(jj==2) omega=-omega ! from GAMESS implementation, this is what seems to be! + if(jj==2) omega=-omega ! from GAMESS implementation, this is what seems to be! call cpu_time(start_time) allocate(inv_resp(nci*(nci+1)/2)) inv_resp=apb-omega**2.0*inv_amb - + uplo='U' allocate(ipiv(1:nci),work(1:nci)) call ssptrf(uplo,nci,inv_resp,ipiv,info) - call ssptri(uplo,nci,inv_resp,ipiv,work,info) + call ssptri(uplo,nci,inv_resp,ipiv,work,info) deallocate(ipiv,work) - - - + + + !$omp parallel private(i,j,jk,io,iv,idum1,idum2,ij) !$omp& reduction (+:XpY) !$omp do @@ -314,9 +314,9 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) io=iconf(j,1) iv=iconf(j,2) idum1=max(io,iv) - idum2=min(io,iv) + idum2=min(io,iv) ij=idum2+idum1*(idum1-1)/2 - alpha_xx=alpha_xx-(xl(ij)*XpY(j,ii,jj,1))*2.0 + alpha_xx=alpha_xx-(xl(ij)*XpY(j,ii,jj,1))*2.0 alpha_xy=alpha_xy-(xl(ij)*XpY(j,ii,jj,2))*2.0 alpha_xz=alpha_xz-(xl(ij)*XpY(j,ii,jj,3))*2.0 alpha_yy=alpha_yy-(yl(ij)*XpY(j,ii,jj,2))*2.0 @@ -337,11 +337,11 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) write(*,111) 'Mean of alpha',(alpha_xx+alpha_yy+alpha_zz)/3.0 write(*,*) endif - - - + + + c frequency-dependent first hyperpolarizability (beta) SHG only - + ! extract X and Y from XpY ! (X-Y)=omega*(A-B)^-1 (X+Y) ! X=((X+Y)+(X-Y))/2 @@ -370,10 +370,10 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) enddo enddo !$omp end do -!$omp end parallel +!$omp end parallel deallocate(inv_resp) - + call cpu_time(end_time) print '("alpha Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 @@ -381,7 +381,7 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) c beta beta = A - B + C (here C=0 since no derivative of g^XC) call cpu_time(start_time) - + if(ii==1)then ! ! Genarating a list of indexes used in A and B formula to save a great bunch of time @@ -407,7 +407,7 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) !Do i=1,counter_A !write(*,*)i,A_list(i,1:3) !enddo - + counter_B=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_B) @@ -455,7 +455,7 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) . A_list,B_list,counter_A,counter_B, . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) - beta(3,3,3)=A-B + beta(3,3,3)=A-B !Off-diagonal !ijj call beta_resp_fast(1,2,2,X,Y, @@ -511,7 +511,7 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) beta(3,2,1)=beta(1,2,3) beta(2,3,1)=beta(1,2,3) beta(2,1,3)=beta(1,2,3) - + else ! ijk=ikj because SHG (2 identical frequencies) !Diagonal call beta_resp_fast(1,1,1,X,Y, @@ -528,7 +528,7 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) . A_list,B_list,counter_A,counter_B, . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) - beta(3,3,3)=A-B + beta(3,3,3)=A-B !Off-diagonal !ijj call beta_resp_fast(1,2,2,X,Y, @@ -542,7 +542,7 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) beta(2,1,2)=A-B beta(2,2,1)=beta(2,1,2) - + call beta_resp_fast(1,3,3,X,Y, . A_list,B_list,counter_A,counter_B, . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) @@ -617,13 +617,13 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) beta(2,3,1)=A-B beta(2,1,3)=beta(2,3,1) - + endif call cpu_time(end_time) print '("beta Time = ",f12.2," minutes.")' - . ,(end_time-start_time)/60.0 - + . ,(end_time-start_time)/60.0 + betaVec_X=0.0 betaVec_Y=0.0 betaVec_Z=0.0 @@ -635,8 +635,8 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) betaVec_X=betaVec_X/5.0 betaVec_Y=betaVec_Y/5.0 betaVec_Z=betaVec_Z/5.0 - - + + write(*,*) '_____________________________________________ .________' write(*,*) @@ -667,167 +667,167 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) .________' write(555,*)freq(ii)*27.21139,beta_HRS - -! beta(:,:,:)=0.0 + +! beta(:,:,:)=0.0 ! ! Let's save time applying symmetry conditions ! if(ii==1)then !static case Kleinman condition applied ijk = kij = jki ! ! only 10 tensor components computed ! !Diagonal ! call beta_resp(1,1,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,1,1)=A-B ! call beta_resp(2,2,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,2,2)=A-B ! call beta_resp(3,3,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! -! beta(3,3,3)=A-B +! +! beta(3,3,3)=A-B ! !Off-diagonal ! !ijj ! call beta_resp(1,2,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,2,2)=A-B ! beta(2,1,2)=beta(1,2,2) ! beta(2,2,1)=beta(1,2,2) ! call beta_resp(1,3,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,3,3)=A-B ! beta(3,1,3)=beta(1,3,3) ! beta(3,3,1)=beta(1,3,3) ! call beta_resp(2,1,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,1,1)=A-B ! beta(1,2,1)=beta(2,1,1) ! beta(1,1,2)=beta(2,1,1) ! call beta_resp(2,3,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,3,3)=A-B ! beta(3,2,3)=beta(2,3,3) ! beta(3,3,2)=beta(2,3,3) ! call beta_resp(3,2,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(3,2,2)=A-B ! beta(2,3,2)=beta(3,2,2) ! beta(2,2,3)=beta(3,2,2) ! call beta_resp(3,1,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(3,1,1)=A-B ! beta(1,3,1)=beta(3,1,1) ! beta(1,1,3)=beta(3,1,1) ! !ijk ! call beta_resp(1,2,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,2,3)=A-B ! beta(1,3,2)=beta(1,2,3) ! beta(3,1,2)=beta(1,2,3) ! beta(3,2,1)=beta(1,2,3) ! beta(2,3,1)=beta(1,2,3) ! beta(2,1,3)=beta(1,2,3) -! +! ! else ! ijk=ikj because SHG (2 identical frequencies) ! !Diagonal ! call beta_resp(1,1,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,1,1)=A-B ! call beta_resp(2,2,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,2,2)=A-B ! call beta_resp(3,3,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! -! beta(3,3,3)=A-B +! +! beta(3,3,3)=A-B ! !Off-diagonal ! !ijj ! call beta_resp(1,2,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,2,2)=A-B ! call beta_resp(2,1,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,1,2)=A-B ! beta(2,2,1)=beta(2,1,2) -! +! ! call beta_resp(1,3,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,3,3)=A-B ! call beta_resp(3,1,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(3,1,3)=A-B ! beta(3,3,1)=beta(3,1,3) ! call beta_resp(2,1,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,1,1)=A-B ! call beta_resp(1,2,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,2,1)=A-B ! beta(1,1,2)=beta(1,2,1) ! call beta_resp(2,3,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,3,3)=A-B ! call beta_resp(3,2,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(3,2,3)=A-B ! beta(3,3,2)=beta(3,2,3) ! call beta_resp(3,2,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(3,2,2)=A-B ! call beta_resp(2,3,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,3,2)=A-B ! beta(2,2,3)=beta(2,3,2) ! call beta_resp(3,1,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(3,1,1)=A-B ! call beta_resp(1,3,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,3,1)=A-B ! beta(1,1,3)=beta(1,3,1) ! !ijk ! call beta_resp(1,2,3,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(1,2,3)=A-B ! beta(1,3,2)=beta(1,2,3) ! call beta_resp(3,1,2,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(3,1,2)=A-B ! beta(3,2,1)=beta(3,1,2) ! call beta_resp(2,3,1,X,Y, ! . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) -! +! ! beta(2,3,1)=A-B ! beta(2,1,3)=beta(2,3,1) -! +! ! endif -! +! ! call cpu_time(end_time) ! print '("beta Time = ",f12.2," minutes.")' -! . ,(end_time-start_time)/60.0 -! +! . ,(end_time-start_time)/60.0 +! ! betaVec_X=0.0 ! betaVec_Y=0.0 ! betaVec_Z=0.0 @@ -839,8 +839,8 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) ! betaVec_X=betaVec_X/5.0 ! betaVec_Y=betaVec_Y/5.0 ! betaVec_Z=betaVec_Z/5.0 -! -! +! +! ! write(*,*) '_____________________________________________ ! .________' ! write(*,*) @@ -859,15 +859,15 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) ! write(*,4444) 'beta_HRS',beta_HRS,' DR', DR ! write(*,*) '_____________________________________________ ! .________' - + enddo deallocate(XpY) deallocate(XmY) deallocate(X,Y) deallocate(inv_amb) -111 format(A15,F20.6) +111 format(A15,F20.6) 1111 format(A22,2A20) -2222 format(A2,3F20.6) +2222 format(A2,3F20.6) 3333 format(A16,F7.1,A1,F7.1,A1) 4444 format(A12,F20.3,A12,F20.3) 5555 format(A31,F7.1,A1,F7.1,A1,F7.1,A1) @@ -877,14 +877,14 @@ SUBROUTINE lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv) . /,2X,'rho=|Bj=3|/|Bj=1| :', F15.4,/) end subroutine lresp1 - + SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, . A_list,B_list,counter_A,counter_B, . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) use commonresp use omp_lib implicit none - + integer ::xx,yy,zz integer ::ix,iy,iz,wx,wy,wz,no,nv,nci,ii,maxconf,moci integer ::iconf(maxconf,2) @@ -893,12 +893,12 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, real*8 ::Y(nci,num_freq+1,2,3) real*8 ::A,B real*8 ::A1,A2,A3,A4,A5,A6 - real*8 ::B1,B2,B3,B4,B5,B6 + real*8 ::B1,B2,B3,B4,B5,B6 integer ::counter_A,counter_B integer :: A_list(1:counter_A,1:3) integer :: B_list(1:counter_B,1:3) - + A=0.0 B=0.0 A1=0.0 @@ -934,7 +934,7 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, call A_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,A2,A_list,counter_A) -c A3 iy iz ix +c A3 iy iz ix xx=iy wx=1 yy=iz @@ -943,8 +943,8 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, wz=2 call A_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,A3,A_list,counter_A) - -c A4 iy ix iz + +c A4 iy ix iz xx=iy wx=1 yy=ix @@ -952,9 +952,9 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, zz=iz wz=1 call A_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, - . mu,nci,moci,ii,A4,A_list,counter_A) - -c A5 iz ix iy + . mu,nci,moci,ii,A4,A_list,counter_A) + +c A5 iz ix iy xx=iz wx=1 yy=ix @@ -963,25 +963,25 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, wz=1 call A_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,A5,A_list,counter_A) - -c A6 iz iy ix + +c A6 iz iy ix xx=iz wx=1 yy=iy wy=1 zz=ix - wz=2 + wz=2 call A_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,A6,A_list,counter_A) - - A=A1+A2+A3+A4+A5+A6 - + + A=A1+A2+A3+A4+A5+A6 + !write(*,*) '______________' !write(*,*) 'A',ix,iy,iz !write(*,*) A1,A2,A3 !write(*,*) A4,A5,A6 !write(*,*) A - + c B1 ix iy iz xx=ix wx=2 @@ -1001,8 +1001,8 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, call B_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,B2,B_list,counter_B) - -c B3 iy iz ix + +c B3 iy iz ix xx=iy wx=1 yy=iz @@ -1012,8 +1012,8 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, call B_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,B3,B_list,counter_B) - -c B4 iy ix iz + +c B4 iy ix iz xx=iy wx=1 yy=ix @@ -1023,8 +1023,8 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, call B_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,B4,B_list,counter_B) - -c B5 iz ix iy + +c B5 iz ix iy xx=iz wx=1 yy=ix @@ -1034,7 +1034,7 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, call B_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,B5,B_list,counter_B) - + c B6 iz iy ix xx=iz wx=1 @@ -1044,18 +1044,18 @@ SUBROUTINE beta_resp_fast(ix,iy,iz,X,Y, wz=2 call B_beta_fast(xx,yy,zz,wx,wy,wz,X,Y, . mu,nci,moci,ii,B6,B_list,counter_B) - + B=B1+B2+B3+B4+B5+B6 - + !write(*,*) 'B',ix,iy,iz !write(*,*) B1,B2,B3 !write(*,*) B4,B5,B6 !write(*,*) B !write(*,*) '______________' - + end subroutine beta_resp_fast - + Subroutine List_A(maxconf,no,nci,iconf,A_list,counter_A) use commonresp use omp_lib @@ -1067,7 +1067,7 @@ Subroutine List_A(maxconf,no,nci,iconf,A_list,counter_A) integer ::A_list(1:counter_A,1:3) integer ::i,j,ij,kk,io1,io2,idum1,idum2 integer ::jwrk - + counter=1 !$omp parallel private(j,i,kk,jwrk,io1,io2,idum1,idum2,ij) !$omp do ordered @@ -1094,9 +1094,9 @@ Subroutine List_A(maxconf,no,nci,iconf,A_list,counter_A) !$OMP END ORDERED enddo !$omp end do -!$omp end parallel +!$omp end parallel end subroutine List_A - + Subroutine List_B(maxconf,no,nv,nci,iconf,B_list,counter_B) use commonresp use omp_lib @@ -1133,10 +1133,10 @@ Subroutine List_B(maxconf,no,nv,nci,iconf,B_list,counter_B) !$OMP END ORDERED enddo !$omp end do -!$omp end parallel +!$omp end parallel end subroutine List_B - - + + Subroutine A_beta_fast(ix,iy,iz,wx,wy,wz,X,Y, . mu,nci,moci,ifreq,A,A_list,counter_A) use commonresp @@ -1149,12 +1149,12 @@ Subroutine A_beta_fast(ix,iy,iz,wx,wy,wz,X,Y, real*8 ::X(nci,num_freq+1,2,3) real*8 ::Y(nci,num_freq+1,2,3) integer ::i,ii - + integer ::counter_A integer ::A_list(1:counter_A,1:3) - + ii=ifreq - A=0.0 + A=0.0 !$omp parallel private(i) !$omp& reduction(+:A) !$omp do @@ -1163,9 +1163,9 @@ Subroutine A_beta_fast(ix,iy,iz,wx,wy,wz,X,Y, . *Y(A_list(i,3),ii,wz,iz) enddo !$omp end do -!$omp end parallel +!$omp end parallel end subroutine A_beta_fast - + Subroutine B_beta_fast(ix,iy,iz,wx,wy,wz,X,Y, . mu,nci,moci,ifreq,B,B_list,counter_B) use commonresp @@ -1178,12 +1178,12 @@ Subroutine B_beta_fast(ix,iy,iz,wx,wy,wz,X,Y, real*8 ::X(nci,num_freq+1,2,3) real*8 ::Y(nci,num_freq+1,2,3) integer ::i,ii - + integer ::counter_B integer ::B_list(1:counter_B,1:3) - + ii=ifreq - B=0.0 + B=0.0 !$omp parallel private(i) !$omp& reduction(+:B) !$omp do @@ -1193,14 +1193,14 @@ Subroutine B_beta_fast(ix,iy,iz,wx,wy,wz,X,Y, enddo !$omp end do !$omp end parallel - end subroutine B_beta_fast - + end subroutine B_beta_fast + SUBROUTINE beta_resp(ix,iy,iz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) use commonresp use omp_lib implicit none - + integer ::xx,yy,zz integer ::ix,iy,iz,wx,wy,wz,no,nv,nci,ii,maxconf,moci integer ::iconf(maxconf,2) @@ -1210,7 +1210,7 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, real*8 ::A,B real*8 ::A1,A2,A3,A4,A5,A6 real*8 ::B1,B2,B3,B4,B5,B6 - + A=0.0 B=0.0 A1=0.0 @@ -1235,7 +1235,7 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, wz=1 call A_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,A1) - + c A2 ix iz iy xx=ix wx=2 @@ -1246,7 +1246,7 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, call A_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,A2) -c A3 iy iz ix +c A3 iy iz ix xx=iy wx=1 yy=iz @@ -1255,7 +1255,7 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, wz=2 call A_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,A3) -c A4 iy ix iz +c A4 iy ix iz xx=iy wx=1 yy=ix @@ -1263,8 +1263,8 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, zz=iz wz=1 call A_beta1(xx,yy,zz,wx,wy,wz,X,Y, - . mu,maxconf,no,nv,nci,moci,ii,iconf,A4) -c A5 iz ix iy + . mu,maxconf,no,nv,nci,moci,ii,iconf,A4) +c A5 iz ix iy xx=iz wx=1 yy=ix @@ -1273,24 +1273,24 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, wz=1 call A_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,A5) -c A6 iz iy ix +c A6 iz iy ix xx=iz wx=1 yy=iy wy=1 zz=ix - wz=2 + wz=2 call A_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,A6) - - A=A1+A2+A3+A4+A5+A6 - + + A=A1+A2+A3+A4+A5+A6 + !write(*,*) '______________' !write(*,*) 'A',ix,iy,iz !write(*,*) A1,A2,A3 !write(*,*) A4,A5,A6 !write(*,*) A - + c B1 ix iy iz xx=ix wx=2 @@ -1300,7 +1300,7 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, wz=1 call B_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,B1) - + c B2 ix iz iy xx=ix wx=2 @@ -1311,8 +1311,8 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, call B_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,B2) - -c B3 iy iz ix + +c B3 iy iz ix xx=iy wx=1 yy=iz @@ -1322,8 +1322,8 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, call B_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,B3) - -c B4 iy ix iz + +c B4 iy ix iz xx=iy wx=1 yy=ix @@ -1333,8 +1333,8 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, call B_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,B4) - -c B5 iz ix iy + +c B5 iz ix iy xx=iz wx=1 yy=ix @@ -1344,7 +1344,7 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, call B_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,B5) - + c B6 iz iy ix xx=iz wx=1 @@ -1354,18 +1354,18 @@ SUBROUTINE beta_resp(ix,iy,iz,X,Y, wz=2 call B_beta1(xx,yy,zz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ii,iconf,B6) - + B=B1+B2+B3+B4+B5+B6 - + !write(*,*) 'B',ix,iy,iz !write(*,*) B1,B2,B3 !write(*,*) B4,B5,B6 !write(*,*) B !write(*,*) '______________' - + end subroutine beta_resp - + Subroutine A_beta1(ix,iy,iz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ifreq,iconf,A) use commonresp @@ -1382,9 +1382,9 @@ Subroutine A_beta1(ix,iy,iz,wx,wy,wz,X,Y, integer*8 ::lin8 integer ::jwrk logical ::check - + ii=ifreq - A=0.0 + A=0.0 !$omp parallel private(j,i,kk,jwrk,io1,io2,idum1,idum2,ij) !$omp& reduction(+:A) !$omp do @@ -1408,9 +1408,9 @@ Subroutine A_beta1(ix,iy,iz,wx,wy,wz,X,Y, enddo enddo !$omp end do -!$omp end parallel +!$omp end parallel end subroutine A_beta1 - + Subroutine B_beta1(ix,iy,iz,wx,wy,wz,X,Y, . mu,maxconf,no,nv,nci,moci,ifreq,iconf,B) use commonresp @@ -1427,9 +1427,9 @@ Subroutine B_beta1(ix,iy,iz,wx,wy,wz,X,Y, integer*8 ::lin8 integer ::jwrk logical ::check - + ii=ifreq - B=0.0 + B=0.0 !$omp parallel private(j,i,kk,jwrk,iv1,iv2,idum1,idum2,ab) !$omp& reduction(+:B) !$omp do @@ -1455,16 +1455,16 @@ Subroutine B_beta1(ix,iy,iz,wx,wy,wz,X,Y, !$omp end do !$omp end parallel end subroutine B_beta1 - + SUBROUTINE HRS(beta,beta2_ZZZ, beta2_XZZ, beta_HRS, DR) implicit none INTEGER::i,j,k REAL*8::beta(1:3,1:3,1:3) - REAL*8:: beta2_ZZZ, beta2_XZZ, beta2_iii, beta2_iij, beta2_jii, - . beta2_ijk, beta_iii_beta_ijj, beta_jii_beta_iij, - . beta_iii_beta_jji,beta_iij_beta_jkk, beta_jii_beta_jkk, - . beta_iij_beta_kkj, beta_ijk_beta_jik, beta2_ijj, - . beta_iij_beta_jii, beta_ijj_beta_ikk, beta_iik_beta_jjk, + REAL*8:: beta2_ZZZ, beta2_XZZ, beta2_iii, beta2_iij, beta2_jii, + . beta2_ijk, beta_iii_beta_ijj, beta_jii_beta_iij, + . beta_iii_beta_jji,beta_iij_beta_jkk, beta_jii_beta_jkk, + . beta_iij_beta_kkj, beta_ijk_beta_jik, beta2_ijj, + . beta_iij_beta_jii, beta_ijj_beta_ikk, beta_iik_beta_jjk, . DR, beta_HRS ! Calcul du beta_HRS @@ -1612,25 +1612,25 @@ SUBROUTINE HRS(beta,beta2_ZZZ, beta2_XZZ, beta_HRS, DR) enddo - beta2_ZZZ=1./7.*beta2_iii + 4./35.*beta2_iij + 1./35.*beta2_jii - . + 2./105.*beta2_ijk + 2./35.*beta_iii_beta_ijj - . + 4./35.*beta_jii_beta_iij + 4./35.*beta_iii_beta_jji - . + 4./105.*beta_iij_beta_jkk + 1./105.*beta_jii_beta_jkk + beta2_ZZZ=1./7.*beta2_iii + 4./35.*beta2_iij + 1./35.*beta2_jii + . + 2./105.*beta2_ijk + 2./35.*beta_iii_beta_ijj + . + 4./35.*beta_jii_beta_iij + 4./35.*beta_iii_beta_jji + . + 4./105.*beta_iij_beta_jkk + 1./105.*beta_jii_beta_jkk . + 4./105.*beta_iij_beta_kkj + 4./105.*beta_ijk_beta_jik - beta2_XZZ=1./35.*beta2_iii + 4./105.*beta_iii_beta_ijj - . - 2./35.*beta_iii_beta_jji + 8./105.*beta2_iij - . - 2./105.*beta_iij_beta_jkk + 2./35.*beta2_ijk - . - 2./105.*beta_ijk_beta_jik + 3./35.*beta2_ijj - . - 2./35.*beta_iij_beta_jii + 1./35.*beta_ijj_beta_ikk + beta2_XZZ=1./35.*beta2_iii + 4./105.*beta_iii_beta_ijj + . - 2./35.*beta_iii_beta_jji + 8./105.*beta2_iij + . - 2./105.*beta_iij_beta_jkk + 2./35.*beta2_ijk + . - 2./105.*beta_ijk_beta_jik + 3./35.*beta2_ijj + . - 2./35.*beta_iij_beta_jii + 1./35.*beta_ijj_beta_ikk . - 2./105.*beta_iik_beta_jjk - + DR=beta2_ZZZ/beta2_XZZ beta_HRS=sqrt(beta2_ZZZ+beta2_XZZ) - end subroutine HRS - + end subroutine HRS + subroutine PrintBeta(beta,unit) implicit none c Arguments @@ -1643,7 +1643,7 @@ subroutine PrintBeta(beta,unit) write(*,9000) (FIELDDIR(i), i=1,3) do i = 1,3 do j = 1,3 - write(*,9001) FIELDDIR(i),FIELDDIR(j), + write(*,9001) FIELDDIR(i),FIELDDIR(j), & (beta(i,j,k), k=1,3) write(unit,9002)(beta(i,j,k), k=1,3) end do @@ -1652,15 +1652,15 @@ subroutine PrintBeta(beta,unit) 9001 format(3X,A1,A1,".",3(2X,F15.6)) 9002 format(3(2X,F15.6)) end - - - + + + SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, - .no,nv,eci,Xci,Yci,nroot) + .no,nv,eci,Xci,Yci,nroot) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot integer ::maxconf,moci,no,nv @@ -1668,13 +1668,13 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, integer, allocatable :: A_list(:,:) integer ::counter_A integer, allocatable :: B_list(:,:) - integer ::counter_B + integer ::counter_B integer*8 ::lin8 - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) real*4 ::omega real*4 ::Xci(nci,nroot), Yci(nci,nroot),eci(nci) @@ -1690,24 +1690,24 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, integer ::info integer, allocatable ::ipiv(:) real*4, allocatable ::work (:) - + integer ::ix,iy,iz real*8 ::sigma(3,3),A,B,sigma_f,sigma_g,sigma_h - + real*8 ::alpha_xx,alpha_xy,alpha_xz real*8 ::alpha_yy,alpha_yz real*8 ::alpha_zz - + real*4 ::start_time,end_time - + open(unit=60,file='2PA-abs',status='replace') - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) mu(:,3)=zl(:) - + write(*,*) write(*,*)'==================================================== .==================' @@ -1716,15 +1716,15 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, write(*,*)'==================================================== .==================' write(*,*) - + allocate(inv_amb(nci*(nci+1)/2)) inv_amb=amb uplo='U' allocate(ipiv(1:nci),work(1:nci)) call ssptrf(uplo,nci,inv_amb,ipiv,info) - call ssptri(uplo,nci,inv_amb,ipiv,work,info) + call ssptri(uplo,nci,inv_amb,ipiv,work,info) deallocate(ipiv,work) - + allocate( XpY(nci,3)) allocate( XmY(nci,3)) allocate( X(nci,3), @@ -1739,16 +1739,16 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, call cpu_time(start_time) allocate(inv_resp(nci*(nci+1)/2)) inv_resp=apb-omega**2.0*inv_amb - + uplo='U' allocate(ipiv(1:nci),work(1:nci)) call ssptrf(uplo,nci,inv_resp,ipiv,info) - call ssptri(uplo,nci,inv_resp,ipiv,work,info) + call ssptri(uplo,nci,inv_resp,ipiv,work,info) deallocate(ipiv,work) - - - + + + !$omp parallel private(i,j,jk,io,iv,idum1,idum2,ij) !$omp& reduction (+:XpY) !$omp do @@ -1784,9 +1784,9 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! io=iconf(j,1) ! iv=iconf(j,2) ! idum1=max(io,iv) -! idum2=min(io,iv) +! idum2=min(io,iv) ! ij=idum2+idum1*(idum1-1)/2 -! alpha_xx=alpha_xx-(xl(ij)*XpY(j,1))*2.0 +! alpha_xx=alpha_xx-(xl(ij)*XpY(j,1))*2.0 ! alpha_xy=alpha_xy-(xl(ij)*XpY(j,2))*2.0 ! alpha_xz=alpha_xz-(xl(ij)*XpY(j,3))*2.0 ! alpha_yy=alpha_yy-(yl(ij)*XpY(j,2))*2.0 @@ -1809,7 +1809,7 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! write(*,111) 'Mean of alpha',(alpha_xx+alpha_yy+alpha_zz)/3.0 ! write(*,*) - + ! extract X and Y from XpY ! (X-Y)=omega*(A-B)^-1 (X+Y) ! X=((X+Y)+(X-Y))/2 @@ -1838,14 +1838,14 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, enddo enddo !$omp end do -!$omp end parallel +!$omp end parallel deallocate(inv_resp) - + call cpu_time(end_time) print '("alpha Time = ",f12.2," minutes.")' - . ,(end_time-start_time)/60.0 - + . ,(end_time-start_time)/60.0 + if(ii==1)then ! ! Genarating a list of indexes used in A and B formula to save a great bunch of time @@ -1872,7 +1872,7 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, !Do i=1,counter_A !write(*,*)i,A_list(i,1:3) !enddo - + counter_B=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_B) @@ -1899,7 +1899,7 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, . ,(end_time-start_time)/60.0 endif -c sigma sigma = -A + B +c sigma sigma = -A + B call cpu_time(start_time) ! ! Fast version @@ -1916,7 +1916,7 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, endif enddo enddo - + sigma_f=0.0 sigma_g=0.0 sigma_h=0.0 @@ -1930,12 +1930,12 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, sigma_f=sigma_f/30.0 sigma_g=sigma_g/30.0 sigma_h=sigma_h/30.0 - - - + + + call cpu_time(end_time) print '("2PA Time = ",f12.2," minutes.")' - . ,(end_time-start_time)/60.0 + . ,(end_time-start_time)/60.0 write(*,*) write(*,3333)'Delta (',eci(ii)*27.21139,')' write(*,*) @@ -1953,9 +1953,9 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, . +3.0*sigma_h write(*,4444)'rho = //*(_|_)**-1 =',(2.0*sigma_f+2.0*sigma_g . +2.0*sigma_h)/(-1.0*sigma_f+4.0*sigma_g-1.0*sigma_h) - + ! call cpu_time(start_time) -! ! +! ! ! ! regular version ! ! ! sigma(:,:)=0.0 @@ -1969,7 +1969,7 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! endif ! enddo ! enddo -! +! ! sigma_f=0.0 ! sigma_g=0.0 ! sigma_h=0.0 @@ -1983,12 +1983,12 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! sigma_f=sigma_f/30.0 ! sigma_g=sigma_g/30.0 ! sigma_h=sigma_h/30.0 -! -! -! +! +! +! ! call cpu_time(end_time) ! print '("2PA Time = ",f12.2," minutes.")' -! . ,(end_time-start_time)/60.0 +! . ,(end_time-start_time)/60.0 ! write(*,*) ! write(*,3333)'Delta (',eci(ii)*27.21139,')' ! write(*,*) @@ -1998,7 +1998,7 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! write(*,2222)'z',sigma(3,1),sigma(3,2),sigma(3,3) ! write(*,*) ! write(*,5555)'F =',sigma_f,' G =',sigma_G,' H =',sigma_H -! !write(*,*)'like in gamess, only F and G, and wrong +! !write(*,*)'like in gamess, only F and G, and wrong ! !.definition for circ. and rho' ! !write(*,*)'Sigma_2PA_// =',2.0*sigma_f+4.0*sigma_g ! !write(*,*)'Sigma_2PA_circ =',-2.0*sigma_f+6.0*sigma_g @@ -2014,7 +2014,7 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . +3.0*sigma_h ! write(*,4444)'rho = //*(_|_)**-1 =',(2.0*sigma_f+2.0*sigma_g ! . +2.0*sigma_h)/(-1.0*sigma_f+4.0*sigma_g-1.0*sigma_h) - + write(60,6666)ii,eci(ii)*27.21139,2.0*sigma_f+2.0*sigma_g .+2.0*sigma_h,-1.0*sigma_f+4.0*sigma_g-1.0*sigma_h,-2.0*sigma_f .+3.0*sigma_g+3.0*sigma_h,(2.0*sigma_f+2.0*sigma_g @@ -2033,16 +2033,16 @@ SUBROUTINE lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, write(*,*)'==================================================== .==================' write(*,*) -111 format(A15,F20.6) +111 format(A15,F20.6) 1111 format(A22,2A20) -2222 format(A2,3F20.6) +2222 format(A2,3F20.6) 3334 format(A16,F7.3,A1,F7.3,A1) 3333 format(A16,F7.3,A1) 4444 format(A20,F20.3) 5555 format(A3,F20.3,A4,F20.3,A4,F20.3) 6666 format(I3,F7.3,4F20.3) end subroutine lresp_2PA - + SUBROUTINE TPA_resp_fast(ix,iy,X,Y,Xci,Yci,nroot, . A_list,B_list,counter_A,counter_B, . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) @@ -2050,7 +2050,7 @@ SUBROUTINE TPA_resp_fast(ix,iy,X,Y,Xci,Yci,nroot, use commonresp use omp_lib implicit none - + integer ::xx,yy,nroot integer ::ix,iy,no,nv,nci,ii,maxconf,moci integer ::iconf(maxconf,2) @@ -2064,7 +2064,7 @@ SUBROUTINE TPA_resp_fast(ix,iy,X,Y,Xci,Yci,nroot, integer ::counter_A,counter_B integer :: A_list(1:counter_A,1:3) integer :: B_list(1:counter_B,1:3) - + A=0.0 B=0.0 A1=0.0 @@ -2101,7 +2101,7 @@ SUBROUTINE TPA_resp_fast(ix,iy,X,Y,Xci,Yci,nroot, . A_list,counter_A, . mu,nci,moci,ii,A4) A=A1+A2+A3+A4 - + c B1 ix iy n xx=ix yy=iy @@ -2127,9 +2127,9 @@ SUBROUTINE TPA_resp_fast(ix,iy,X,Y,Xci,Yci,nroot, . B_list,counter_B, . mu,nci,moci,ii,B4) B=B1+B2+B3+B4 - + end subroutine TPA_resp_fast - + Subroutine A_2PA_1_fast(ix,iy,X,Y,Xci,Yci,nroot, . A_list,counter_A, . mu,nci,moci,ifreq,A) @@ -2145,12 +2145,12 @@ Subroutine A_2PA_1_fast(ix,iy,X,Y,Xci,Yci,nroot, real*8 ::Y(nci,3) real*4 ::Xci(nci,nroot), Yci(nci,nroot) integer ::i,ii - + integer ::counter_A integer ::A_list(1:counter_A,1:3) - + ii=ifreq - A=0.0 + A=0.0 !$omp parallel private(i) !$omp& reduction(+:A) !$omp do @@ -2159,9 +2159,9 @@ Subroutine A_2PA_1_fast(ix,iy,X,Y,Xci,Yci,nroot, . *dble(Yci(A_list(i,3),ii)) enddo !$omp end do -!$omp end parallel +!$omp end parallel end subroutine A_2PA_1_fast - + Subroutine A_2PA_2_fast(ix,iy,X,Y,Xci,Yci,nroot, . A_list,counter_A, . mu,nci,moci,ifreq,A) @@ -2177,12 +2177,12 @@ Subroutine A_2PA_2_fast(ix,iy,X,Y,Xci,Yci,nroot, real*8 ::Y(nci,3) real*4 ::Xci(nci,nroot), Yci(nci,nroot) integer ::i,ii - + integer ::counter_A integer ::A_list(1:counter_A,1:3) - + ii=ifreq - A=0.0 + A=0.0 !$omp parallel private(i) !$omp& reduction(+:A) !$omp do @@ -2191,9 +2191,9 @@ Subroutine A_2PA_2_fast(ix,iy,X,Y,Xci,Yci,nroot, . *Y(A_list(i,3),iy) enddo !$omp end do -!$omp end parallel +!$omp end parallel end subroutine A_2PA_2_fast - + Subroutine B_2PA_1_fast(ix,iy,X,Y,Xci,Yci,nroot, . B_list,counter_B, . mu,nci,moci,ifreq,B) @@ -2209,12 +2209,12 @@ Subroutine B_2PA_1_fast(ix,iy,X,Y,Xci,Yci,nroot, real*8 ::Y(nci,3) real*4 ::Xci(nci,nroot), Yci(nci,nroot) integer ::i,ii - + integer ::counter_B integer ::B_list(1:counter_B,1:3) - + ii=ifreq - B=0.0 + B=0.0 !$omp parallel private(i) !$omp& reduction(+:B) !$omp do @@ -2225,7 +2225,7 @@ Subroutine B_2PA_1_fast(ix,iy,X,Y,Xci,Yci,nroot, !$omp end do !$omp end parallel end subroutine B_2PA_1_fast - + Subroutine B_2PA_2_fast(ix,iy,X,Y,Xci,Yci,nroot, . B_list,counter_B, . mu,nci,moci,ifreq,B) @@ -2241,12 +2241,12 @@ Subroutine B_2PA_2_fast(ix,iy,X,Y,Xci,Yci,nroot, real*8 ::Y(nci,3) real*4 ::Xci(nci,nroot), Yci(nci,nroot) integer ::i,ii - + integer ::counter_B integer ::B_list(1:counter_B,1:3) - + ii=ifreq - B=0.0 + B=0.0 !$omp parallel private(i) !$omp& reduction(+:B) !$omp do @@ -2256,16 +2256,16 @@ Subroutine B_2PA_2_fast(ix,iy,X,Y,Xci,Yci,nroot, enddo !$omp end do !$omp end parallel - end subroutine B_2PA_2_fast - - - + end subroutine B_2PA_2_fast + + + SUBROUTINE TPA_resp(ix,iy,X,Y,Xci,Yci,nroot, . mu,maxconf,no,nv,nci,moci,ii,iconf,A,B) use commonresp use omp_lib implicit none - + integer ::xx,yy,nroot integer ::ix,iy,no,nv,nci,ii,maxconf,moci integer ::iconf(maxconf,2) @@ -2276,7 +2276,7 @@ SUBROUTINE TPA_resp(ix,iy,X,Y,Xci,Yci,nroot, real*8 ::A,B real*8 ::A1,A2,A3,A4 real*8 ::B1,B2,B3,B4 - + A=0.0 B=0.0 A1=0.0 @@ -2309,7 +2309,7 @@ SUBROUTINE TPA_resp(ix,iy,X,Y,Xci,Yci,nroot, call A_2PA_2(xx,yy,X,Y,Xci,Yci,nroot, . mu,maxconf,no,nv,nci,moci,ii,iconf,A4) A=A1+A2+A3+A4 - + c B1 ix iy n xx=ix yy=iy @@ -2331,12 +2331,12 @@ SUBROUTINE TPA_resp(ix,iy,X,Y,Xci,Yci,nroot, call B_2PA_2(xx,yy,X,Y,Xci,Yci,nroot, . mu,maxconf,no,nv,nci,moci,ii,iconf,B4) B=B1+B2+B3+B4 - + end subroutine TPA_resp - - - - + + + + Subroutine A_2PA_1(ix,iy,X,Y,Xci,Yci,nroot, . mu,maxconf,no,nv,nci,moci,ifreq,iconf,A) use commonresp @@ -2354,9 +2354,9 @@ Subroutine A_2PA_1(ix,iy,X,Y,Xci,Yci,nroot, integer*8 ::lin8 integer ::jwrk logical ::check - + ii=ifreq - A=0.0 + A=0.0 !$omp parallel private(j,i,kk,jwrk,io1,io2,idum1,idum2,ij) !$omp& reduction(+:A) !$omp do @@ -2380,9 +2380,9 @@ Subroutine A_2PA_1(ix,iy,X,Y,Xci,Yci,nroot, enddo enddo !$omp end do -!$omp end parallel +!$omp end parallel end subroutine A_2PA_1 - + Subroutine A_2PA_2(ix,iy,X,Y,Xci,Yci,nroot, . mu,maxconf,no,nv,nci,moci,ifreq,iconf,A) use commonresp @@ -2400,9 +2400,9 @@ Subroutine A_2PA_2(ix,iy,X,Y,Xci,Yci,nroot, integer*8 ::lin8 integer ::jwrk logical ::check - + ii=ifreq - A=0.0 + A=0.0 !$omp parallel private(j,i,kk,jwrk,io1,io2,idum1,idum2,ij) !$omp& reduction(+:A) !$omp do @@ -2426,9 +2426,9 @@ Subroutine A_2PA_2(ix,iy,X,Y,Xci,Yci,nroot, enddo enddo !$omp end do -!$omp end parallel +!$omp end parallel end subroutine A_2PA_2 - + Subroutine B_2PA_1(ix,iy,X,Y,Xci,Yci,nroot, . mu,maxconf,no,nv,nci,moci,ifreq,iconf,B) use commonresp @@ -2446,9 +2446,9 @@ Subroutine B_2PA_1(ix,iy,X,Y,Xci,Yci,nroot, integer*8 ::lin8 integer ::jwrk logical ::check - + ii=ifreq - B=0.0 + B=0.0 !$omp parallel private(j,i,kk,jwrk,iv1,iv2,idum1,idum2,ab) !$omp& reduction(+:B) !$omp do @@ -2474,7 +2474,7 @@ Subroutine B_2PA_1(ix,iy,X,Y,Xci,Yci,nroot, !$omp end do !$omp end parallel end subroutine B_2PA_1 - + Subroutine B_2PA_2(ix,iy,X,Y,Xci,Yci,nroot, . mu,maxconf,no,nv,nci,moci,ifreq,iconf,B) use commonresp @@ -2492,9 +2492,9 @@ Subroutine B_2PA_2(ix,iy,X,Y,Xci,Yci,nroot, integer*8 ::lin8 integer ::jwrk logical ::check - + ii=ifreq - B=0.0 + B=0.0 !$omp parallel private(j,i,kk,jwrk,iv1,iv2,idum1,idum2,ab) !$omp& reduction(+:B) !$omp do @@ -2520,38 +2520,38 @@ Subroutine B_2PA_2(ix,iy,X,Y,Xci,Yci,nroot, !$omp end do !$omp end parallel end subroutine B_2PA_2 - + subroutine pol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, .maxconf,iconf,ak) use commonresp use omp_lib implicit none - + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot integer ::maxconf,moci integer ::iconf(maxconf,2) integer*8 ::lin8 - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) real*4 ::omega real*4 ::freq(num_freq+1) real*4 ::Xci(nci,nroot), Yci(nci,nroot),eci(nci) - + real*8 ::alpha_xx,alpha_xy,alpha_xz real*8 ::alpha_yy,alpha_yz real*8 ::alpha_zz,ak real*8 ::x,y,z - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) mu(:,3)=zl(:) - + open(unit=101,file='wavelength',form='formatted',status='old') freq(1)=0.0 @@ -2573,7 +2573,7 @@ subroutine pol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, alpha_yy=0.0 alpha_yz=0.0 alpha_zz=0.0 - + Do i=1,nroot x=0.0 y=0.0 @@ -2582,14 +2582,14 @@ subroutine pol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, io=iconf(j,1) iv=iconf(j,2) idum1=max(io,iv) - idum2=min(io,iv) + idum2=min(io,iv) ij=idum2+idum1*(idum1-1)/2 - + x=x+mu(ij,1)*(dble(Xci(j,i))+dble(Yci(j,i))) y=y+mu(ij,2)*(dble(Xci(j,i))+dble(Yci(j,i))) z=z+mu(ij,3)*(dble(Xci(j,i))+dble(Yci(j,i))) enddo - + alpha_xx=alpha_xx-(ak*x*x) . *((1.0/(dble(omega)-dble(eci(i)))) . -(1.0/(dble(omega)+dble(eci(i))))) @@ -2598,7 +2598,7 @@ subroutine pol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, . -(1.0/(dble(omega)+dble(eci(i))))) alpha_zz=alpha_zz-(ak*z*z) . *((1.0/(dble(omega)-dble(eci(i)))) - . -(1.0/(dble(omega)+dble(eci(i))))) + . -(1.0/(dble(omega)+dble(eci(i))))) alpha_xy=alpha_xy-(ak*x*y) . *((1.0/(dble(omega)-dble(eci(i)))) . -(1.0/(dble(omega)+dble(eci(i))))) @@ -2620,22 +2620,22 @@ subroutine pol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, write(*,*) write(*,111) 'Mean of alpha',(alpha_xx+alpha_yy+alpha_zz)/3.0 write(*,*) - - - + + + enddo -111 format(A15,F20.6) +111 format(A15,F20.6) 1111 format(A22,2A20) -2222 format(A2,3F20.6) +2222 format(A2,3F20.6) 3333 format(A16,F7.1,A1,F7.1,A1) end subroutine pol_sos - + subroutine lresp_ESA(nci,iconf,maxconf,xl,yl,zl,moci, - .no,nv,eci,Xci,Yci,nroot,xmolw,thr) + .no,nv,eci,Xci,Yci,nroot,xmolw,thr) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot integer ::maxconf,moci,no,nv @@ -2643,31 +2643,31 @@ subroutine lresp_ESA(nci,iconf,maxconf,xl,yl,zl,moci, integer, allocatable :: A_list(:,:) integer ::counter_A integer, allocatable :: B_list(:,:) - integer ::counter_B + integer ::counter_B integer*8 ::lin8 - + real*8 ::xmolw - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) real*4 ::Xci(nci,nroot), Yci(nci,nroot),eci(nci) real*4 ::delta_e(nroot) real*8 ::thr,umerk(4,nroot) - + real*8 ::mu_s2s(nroot,3),A1,A2,B1,B2,osc_strength(nroot) - + integer ::ix,iy,iz - + real*4 ::start_time,end_time - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) mu(:,3)=zl(:) - + write(*,*) write(*,*)'==================================================== .==================' @@ -2676,7 +2676,7 @@ subroutine lresp_ESA(nci,iconf,maxconf,xl,yl,zl,moci, write(*,*)'==================================================== .==================' write(*,*) - + ! ! Genarating a list of indexes used in A and B formula to save a great bunch of time ! @@ -2702,7 +2702,7 @@ subroutine lresp_ESA(nci,iconf,maxconf,xl,yl,zl,moci, !Do i=1,counter_A !write(*,*)i,A_list(i,1:3) !enddo - + counter_B=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_B) @@ -2777,7 +2777,7 @@ subroutine lresp_ESA(nci,iconf,maxconf,xl,yl,zl,moci, delta_e(ii)=eci(ii)-eci(state2opt) osc_strength(ii)=2.0/3.0*dble(delta_e(ii))* .(mu_s2s(ii,1)**2.0+mu_s2s(ii,2)**2.0+mu_s2s(ii,3)**2.0) - + write(*,11) state2opt, ii, (delta_e(ii))*27.21139, .1.d+7/(delta_e(ii)*2.19474625d+5), .osc_strength(ii) @@ -2795,15 +2795,15 @@ subroutine lresp_ESA(nci,iconf,maxconf,xl,yl,zl,moci, write(*,*) 11 format(i6,i9,f9.3,f9.1,f11.4) 12 format(a1,i3,a17,i3,a1) - + end subroutine lresp_ESA - + subroutine lresp_ESAbis(nci,iconf,maxconf,xl,yl,zl,moci, - .no,nv,eci,Xci,Yci,nroot,mu_s2s) + .no,nv,eci,Xci,Yci,nroot,mu_s2s) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot integer ::maxconf,moci,no,nv @@ -2811,27 +2811,27 @@ subroutine lresp_ESAbis(nci,iconf,maxconf,xl,yl,zl,moci, integer, allocatable :: A_list(:,:) integer ::counter_A integer, allocatable :: B_list(:,:) - integer ::counter_B + integer ::counter_B integer*8 ::lin8 - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) real*4 ::Xci(nci,nroot), Yci(nci,nroot),eci(nci) - + real*8 ::mu_s2s(nroot,nroot,3),A1,A2,B1,B2,osc_strength - + integer ::ix,iy,iz - + real*4 ::start_time,end_time - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) mu(:,3)=zl(:) - + write(*,*) write(*,*)'==================================================== .==================' @@ -2840,7 +2840,7 @@ subroutine lresp_ESAbis(nci,iconf,maxconf,xl,yl,zl,moci, write(*,*)'==================================================== .==================' write(*,*) - + ! ! Genarating a list of indexes used in A and B formula to save a great bunch of time ! @@ -2866,7 +2866,7 @@ subroutine lresp_ESAbis(nci,iconf,maxconf,xl,yl,zl,moci, !Do i=1,counter_A !write(*,*)i,A_list(i,1:3) !enddo - + counter_B=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_B) @@ -2956,44 +2956,44 @@ subroutine lresp_ESAbis(nci,iconf,maxconf,xl,yl,zl,moci, write(*,*) 11 format(i6,i9,f9.3,f11.4) 12 format(a1,i3,a17,i3,a1) - + end subroutine lresp_ESAbis - - + + subroutine hyperpol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, .maxconf,iconf,no,nv) use commonresp use omp_lib implicit none - - ! This subroutine uses dipole moments that are unrelaxed, + + ! This subroutine uses dipole moments that are unrelaxed, ! the result should be different that with response functions - + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot,no,nv integer ::ix,iy,iz,iroot integer ::maxconf,moci integer ::iconf(maxconf,2) integer*8 ::lin8 - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) real*4 ::omega real*4 ::freq(num_freq+1) real*4 ::Xci(nci,nroot), Yci(nci,nroot),eci(nci) - + real*8 ::beta(3,3,3) real*8 ::mu_s(nroot,3) real*8 ::mu_s2s(nroot,nroot,3) - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) mu(:,3)=zl(:) - + open(unit=101,file='wavelength',form='formatted',status='old') freq(1)=0.0 num_freq=1 @@ -3013,9 +3013,9 @@ subroutine hyperpol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, io=iconf(j,1) iv=iconf(j,2) idum1=max(io,iv) - idum2=min(io,iv) + idum2=min(io,iv) ij=idum2+idum1*(idum1-1)/2 - + mu_s(i,1)=mu_s(i,1)+mu(ij,1)*(dble(Xci(j,i))+dble(Yci(j,i))) mu_s(i,2)=mu_s(i,2)+mu(ij,2)*(dble(Xci(j,i))+dble(Yci(j,i))) mu_s(i,3)=mu_s(i,3)+mu(ij,3)*(dble(Xci(j,i))+dble(Yci(j,i))) @@ -3029,7 +3029,7 @@ subroutine hyperpol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, Do iroot=1, num_freq+1 omega=freq(iroot) beta=0.0 - + Do ix=1,3 Do iy=1,3 Do iz=1,3 @@ -3071,38 +3071,38 @@ subroutine hyperpol_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, .________' enddo 5555 format(A31,F7.1,A1,F7.1,A1,F7.1,A1) - end subroutine hyperpol_sos - - + end subroutine hyperpol_sos + + subroutine tpa_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, .maxconf,iconf,no,nv) use commonresp use omp_lib implicit none - - ! This subroutine uses dipole moments that are unrelaxed, + + ! This subroutine uses dipole moments that are unrelaxed, ! the result should be different that with response functions - + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot,no,nv integer ::ix,iy,iz,iroot integer ::maxconf,moci integer ::iconf(maxconf,2) integer*8 ::lin8 - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) real*4 ::omega real*4 ::freq(num_freq+1) real*4 ::Xci(nci,nroot), Yci(nci,nroot),eci(nci) - + real*8 ::sigma(3,3) real*8 ::mu_s(nroot,3) real*8 ::mu_s2s(nroot,nroot,3) - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) @@ -3114,23 +3114,23 @@ subroutine tpa_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, io=iconf(j,1) iv=iconf(j,2) idum1=max(io,iv) - idum2=min(io,iv) + idum2=min(io,iv) ij=idum2+idum1*(idum1-1)/2 - + mu_s(i,1)=mu_s(i,1)+mu(ij,1)*(dble(Xci(j,i))+dble(Yci(j,i))) mu_s(i,2)=mu_s(i,2)+mu(ij,2)*(dble(Xci(j,i))+dble(Yci(j,i))) mu_s(i,3)=mu_s(i,3)+mu(ij,3)*(dble(Xci(j,i))+dble(Yci(j,i))) enddo mu_s(i,:)=mu_s(i,:)*2.0**(1.0/2.0) enddo - + call lresp_ESAbis(nci,iconf,maxconf,xl,yl,zl,moci, .no,nv,eci,Xci,Yci,nroot,mu_s2s) - + Do iroot=1, nroot omega=eci(iroot)/2.0 sigma=0.0 - + Do ix=1,3 Do iy=1,3 Do ii=1,nroot @@ -3151,15 +3151,15 @@ subroutine tpa_sos(nroot,nci,eci,Xci,Yci,xl,yl,zl,moci, enddo 3333 format(A16,F7.3,A1) 1111 format(A22,2A20) -2222 format(A2,3F20.6) - end subroutine tpa_sos - +2222 format(A2,3F20.6) + end subroutine tpa_sos + subroutine lresp_ESA_tda(nci,iconf,maxconf,xl,yl,zl,moci, - .no,nv,eci,Xci,nroot,xmolw,thr) + .no,nv,eci,Xci,nroot,xmolw,thr) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot integer ::maxconf,moci,no,nv @@ -3167,31 +3167,31 @@ subroutine lresp_ESA_tda(nci,iconf,maxconf,xl,yl,zl,moci, integer, allocatable :: A_list(:,:) integer ::counter_A integer, allocatable :: B_list(:,:) - integer ::counter_B + integer ::counter_B integer*8 ::lin8 - + real*8 ::xmolw - + real*8 ::xl(moci*(moci+1)/2) real*8 ::yl(moci*(moci+1)/2) real*8 ::zl(moci*(moci+1)/2) - + real*8 ::mu(moci*(moci+1)/2,3) real*4 ::Xci(nci,nroot),eci(nci) real*4 ::delta_e(nroot) real*8 ::thr,umerk(4,nroot) - + real*8 ::mu_s2s(nroot,3),A1,B1,osc_strength(nroot) - + integer ::ix,iy,iz - + real*4 ::start_time,end_time - + mu=0.0 mu(:,1)=xl(:) mu(:,2)=yl(:) mu(:,3)=zl(:) - + write(*,*) write(*,*)'==================================================== .==================' @@ -3200,7 +3200,7 @@ subroutine lresp_ESA_tda(nci,iconf,maxconf,xl,yl,zl,moci, write(*,*)'==================================================== .==================' write(*,*) - + ! ! Genarating a list of indexes used in A and B formula to save a great bunch of time ! @@ -3226,7 +3226,7 @@ subroutine lresp_ESA_tda(nci,iconf,maxconf,xl,yl,zl,moci, !Do i=1,counter_A !write(*,*)i,A_list(i,1:3) !enddo - + counter_B=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_B) @@ -3295,7 +3295,7 @@ subroutine lresp_ESA_tda(nci,iconf,maxconf,xl,yl,zl,moci, delta_e(ii)=eci(ii)-eci(state2opt) osc_strength(ii)=2.0/3.0*dble(delta_e(ii))* .(mu_s2s(ii,1)**2.0+mu_s2s(ii,2)**2.0+mu_s2s(ii,3)**2.0) - + write(*,11) state2opt, ii, (delta_e(ii))*27.21139, .1.d+7/(delta_e(ii)*2.19474625d+5), .osc_strength(ii) @@ -3313,17 +3313,17 @@ subroutine lresp_ESA_tda(nci,iconf,maxconf,xl,yl,zl,moci, write(*,*) 11 format(i6,i9,f9.3,f9.1,f11.4) 12 format(a1,i3,a17,i3,a1) - - end subroutine lresp_ESA_tda - - + + end subroutine lresp_ESA_tda + + subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, . maxconfb,xla,yla,zla,xlb,ylb,zlb,mocia,mocib, - . noa,nob,nva,nvb,eci,Xci,Yci,nroot,xmolw,thr,ak) + . noa,nob,nva,nvb,eci,Xci,Yci,nroot,xmolw,thr,ak) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::nexa,nexb,maxconfa,maxconfb,mocia,mocib integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot @@ -3332,35 +3332,35 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, integer, allocatable :: A_lista(:,:) integer ::counter_Aa integer, allocatable :: B_lista(:,:) - integer ::counter_Ba + integer ::counter_Ba integer, allocatable :: A_listb(:,:) integer ::counter_Ab integer, allocatable :: B_listb(:,:) - integer ::counter_Bb + integer ::counter_Bb integer*8 ::lin8 - + real*8 ::xmolw,ak - + real*8 ::xla(mocia*(mocia+1)/2) real*8 ::yla(mocia*(mocia+1)/2) real*8 ::zla(mocia*(mocia+1)/2) - + real*8 ::xlb(mocib*(mocib+1)/2) real*8 ::ylb(mocib*(mocib+1)/2) real*8 ::zlb(mocib*(mocib+1)/2) - + real*8 ::mua(mocia*(mocia+1)/2,3) real*8 ::mub(mocib*(mocib+1)/2,3) real*4 ::Xci(nci,nroot), Yci(nci,nroot),eci(nci) real*4 ::delta_e(nroot) real*8 ::thr,umerk(4,nroot) - + real*8 ::mu_s2s(nroot,3),Aa1,Aa2,Ba1,Ba2,osc_strength(nroot) real*8 ::Ab1,Ab2,Bb1,Bb2 integer ::ix,iy,iz - + real*4 ::start_time,end_time - + mua=0.0 mub=0.0 mua(:,1)=xla(:) @@ -3369,7 +3369,7 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, mub(:,1)=xlb(:) mub(:,2)=ylb(:) mub(:,3)=zlb(:) - + write(*,*) write(*,*)'==================================================== .==================' @@ -3378,16 +3378,16 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, write(*,*)'==================================================== .==================' write(*,*) - - - + + + ! ! Genarating a list of indexes used in A and B formula to save a great bunch of time ! - + ! alpha block - - + + call cpu_time(start_time) counter_Aa=0 !$omp parallel private(j,i,kk) @@ -3408,7 +3408,7 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, A_lista=-9999 call List_A(maxconfa,noa,nexa,iconfa,A_lista,counter_Aa) - + counter_Ba=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_Ba) @@ -3427,10 +3427,10 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, allocate(B_lista(1:counter_Ba,1:3)) B_lista=-9999 call List_B(maxconfa,noa,nva,nexa,iconfa,B_lista,counter_Ba) - + ! beta block - - + + call cpu_time(start_time) counter_Ab=0 !$omp parallel private(j,i,kk) @@ -3451,7 +3451,7 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, A_listb=-9999 call List_A(maxconfb,nob,nexb,iconfb,A_listb,counter_Ab) - + counter_Bb=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_Bb) @@ -3470,8 +3470,8 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, allocate(B_listb(1:counter_Bb,1:3)) B_listb=-9999 call List_B(maxconfb,nob,nvb,nexb,iconfb,B_listb,counter_Bb) - - + + call cpu_time(end_time) print '("A & B indexes list Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 @@ -3492,10 +3492,10 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, Ab2=0.0 Bb1=0.0 Bb2=0.0 - + ! alpha - - + + !$omp parallel private(i) !$omp& reduction(+:Aa1,Aa2) !$omp do @@ -3521,10 +3521,10 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, !$omp end do !$omp end parallel - + ! beta - - + + !$omp parallel private(i) !$omp& reduction(+:Ab1,Ab2) !$omp do @@ -3552,7 +3552,7 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, ! divided by two because the beta formula is divided by two in the unrestricted case - + mu_s2s(ii,ix)=-(Aa1+Aa2-Ba1-Ba2+Ab1+Ab2-Bb1-Bb2)/2.0 write(*,*) ix, mu_s2s(ii,ix) enddo @@ -3568,7 +3568,7 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, delta_e(ii)=eci(ii)-eci(state2opt) osc_strength(ii)=2.0/3.0*dble(delta_e(ii))* .(mu_s2s(ii,1)**2.0+mu_s2s(ii,2)**2.0+mu_s2s(ii,3)**2.0) - + write(*,11) state2opt, ii, (delta_e(ii))*27.21139, .1.d+7/(delta_e(ii)*2.19474625d+5), .osc_strength(ii) @@ -3586,16 +3586,16 @@ subroutine ulresp_ESA(nexa,nexb,nci,iconfa,iconfb,maxconfa, write(*,*) 11 format(i6,i9,f9.3,f9.1,f11.4) 12 format(a1,i3,a17,i3,a1) - + end subroutine ulresp_ESA - + subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, . maxconfa,maxconfb,xla,yla,zla,xlb,ylb,zlb,mocia,mocib, - . noa,nob,nva,nvb,eci,Xci,nroot,xmolw,thr,ak) + . noa,nob,nva,nvb,eci,Xci,nroot,xmolw,thr,ak) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::nexa,nexb,maxconfa,maxconfb,mocia,mocib integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot @@ -3604,35 +3604,35 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, integer, allocatable :: A_lista(:,:) integer ::counter_Aa integer, allocatable :: B_lista(:,:) - integer ::counter_Ba + integer ::counter_Ba integer, allocatable :: A_listb(:,:) integer ::counter_Ab integer, allocatable :: B_listb(:,:) - integer ::counter_Bb + integer ::counter_Bb integer*8 ::lin8 - + real*8 ::xmolw,ak - + real*8 ::xla(mocia*(mocia+1)/2) real*8 ::yla(mocia*(mocia+1)/2) real*8 ::zla(mocia*(mocia+1)/2) - + real*8 ::xlb(mocib*(mocib+1)/2) real*8 ::ylb(mocib*(mocib+1)/2) real*8 ::zlb(mocib*(mocib+1)/2) - + real*8 ::mua(mocia*(mocia+1)/2,3) real*8 ::mub(mocib*(mocib+1)/2,3) real*4 ::Xci(nci,nroot),eci(nci) real*4 ::delta_e(nroot) real*8 ::thr,umerk(4,nroot) - + real*8 ::mu_s2s(nroot,3),Aa1,Ba1,osc_strength(nroot) real*8 ::Ab1,Bb1 integer ::ix,iy,iz - + real*4 ::start_time,end_time - + mua=0.0 mub=0.0 mua(:,1)=xla(:) @@ -3641,7 +3641,7 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, mub(:,1)=xlb(:) mub(:,2)=ylb(:) mub(:,3)=zlb(:) - + write(*,*) write(*,*)'==================================================== .==================' @@ -3650,16 +3650,16 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, write(*,*)'==================================================== .==================' write(*,*) - - - + + + ! ! Genarating a list of indexes used in A and B formula to save a great bunch of time ! - + ! alpha block - - + + call cpu_time(start_time) counter_Aa=0 !$omp parallel private(j,i,kk) @@ -3680,7 +3680,7 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, A_lista=-9999 call List_A(maxconfa,noa,nexa,iconfa,A_lista,counter_Aa) - + counter_Ba=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_Ba) @@ -3699,10 +3699,10 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, allocate(B_lista(1:counter_Ba,1:3)) B_lista=-9999 call List_B(maxconfa,noa,nva,nexa,iconfa,B_lista,counter_Ba) - + ! beta block - - + + call cpu_time(start_time) counter_Ab=0 !$omp parallel private(j,i,kk) @@ -3723,7 +3723,7 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, A_listb=-9999 call List_A(maxconfb,nob,nexb,iconfb,A_listb,counter_Ab) - + counter_Bb=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_Bb) @@ -3742,8 +3742,8 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, allocate(B_listb(1:counter_Bb,1:3)) B_listb=-9999 call List_B(maxconfb,nob,nvb,nexb,iconfb,B_listb,counter_Bb) - - + + call cpu_time(end_time) print '("A & B indexes list Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 @@ -3760,10 +3760,10 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, Ba1=0.0 Ab1=0.0 Bb1=0.0 - + ! alpha - - + + !$omp parallel private(i) !$omp& reduction(+:Aa1) !$omp do @@ -3783,10 +3783,10 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, !$omp end do !$omp end parallel - + ! beta - - + + !$omp parallel private(i) !$omp& reduction(+:Ab1) !$omp do @@ -3823,7 +3823,7 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, delta_e(ii)=eci(ii)-eci(state2opt) osc_strength(ii)=2.0/3.0*dble(delta_e(ii))* .(mu_s2s(ii,1)**2.0+mu_s2s(ii,2)**2.0+mu_s2s(ii,3)**2.0) - + write(*,11) state2opt, ii, (delta_e(ii))*27.21139, .1.d+7/(delta_e(ii)*2.19474625d+5), .osc_strength(ii) @@ -3841,16 +3841,16 @@ subroutine ulresp_ESA_tda(nexa,nexb,nci,iconfa,iconfb, write(*,*) 11 format(i6,i9,f9.3,f9.1,f11.4) 12 format(a1,i3,a17,i3,a1) - + end subroutine ulresp_ESA_tda - + subroutine sf_lresp_ESA(nci,iconf,maxconf,xla,yla,zla,mocia, . xlb,ylb,zlb,mocib, - . noa,nva,nob,nvb,eci,Xci,nroot,xmolw,thr) + . noa,nva,nob,nvb,eci,Xci,nroot,xmolw,thr) use commonresp use omp_lib - IMPLICIT NONE - + IMPLICIT NONE + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci integer ::io1,io2,iv1,iv2,iwrk,jwrk,nroot integer ::maxconf,mocia,noa,nva,mocib,nob,nvb @@ -3858,29 +3858,29 @@ subroutine sf_lresp_ESA(nci,iconf,maxconf,xla,yla,zla,mocia, integer, allocatable :: A_list(:,:) integer ::counter_A integer, allocatable :: B_list(:,:) - integer ::counter_B + integer ::counter_B integer*8 ::lin8 - + real*8 ::xmolw - + real*8 ::xla(mocia*(mocia+1)/2) real*8 ::yla(mocia*(mocia+1)/2) real*8 ::zla(mocia*(mocia+1)/2) - + real*8 ::xlb(mocib*(mocib+1)/2) real*8 ::ylb(mocib*(mocib+1)/2) real*8 ::zlb(mocib*(mocib+1)/2) - + real*8 ::mua(mocia*(mocia+1)/2,3) - real*8 ::mub(mocib*(mocib+1)/2,3) + real*8 ::mub(mocib*(mocib+1)/2,3) real*4 ::Xci(nci,nroot),eci(nci) real*4 ::delta_e(nroot) real*8 ::thr,umerk(4,nroot) - + real*8 ::mu_s2s(nroot,3),A1,B1,osc_strength(nroot) - + integer ::ix,iy,iz - + real*4 ::start_time,end_time !first spinflip state @@ -3894,7 +3894,7 @@ subroutine sf_lresp_ESA(nci,iconf,maxconf,xla,yla,zla,mocia, mub(:,1)=xlb(:) mub(:,2)=ylb(:) mub(:,3)=zlb(:) - + write(*,*) write(*,*)'==================================================== .==================' @@ -3904,7 +3904,7 @@ subroutine sf_lresp_ESA(nci,iconf,maxconf,xla,yla,zla,mocia, .==================' write(*,*) write(*,*)'Useful if the first SF state is the ground state' - write(*,*) + write(*,*) ! ! Genarating a list of indexes used in A and B formula to save a great bunch of time ! @@ -3930,7 +3930,7 @@ subroutine sf_lresp_ESA(nci,iconf,maxconf,xla,yla,zla,mocia, !Do i=1,counter_A !write(*,*)i,A_list(i,1:3) !enddo - + counter_B=0 !$omp parallel private(j,i,kk) !$omp& reduction(+:counter_B) @@ -3999,7 +3999,7 @@ subroutine sf_lresp_ESA(nci,iconf,maxconf,xla,yla,zla,mocia, delta_e(ii)=eci(ii)-eci(state2opt) osc_strength(ii)=2.0/3.0*dble(delta_e(ii))* .(mu_s2s(ii,1)**2.0+mu_s2s(ii,2)**2.0+mu_s2s(ii,3)**2.0) - + write(*,11) state2opt, ii, (delta_e(ii))*27.21139, .1.d+7/(delta_e(ii)*2.19474625d+5), .osc_strength(ii) @@ -4017,5 +4017,524 @@ subroutine sf_lresp_ESA(nci,iconf,maxconf,xla,yla,zla,mocia, write(*,*) 11 format(i6,i9,f9.3,f9.1,f11.4) 12 format(a1,i3,a17,i3,a1) - + end subroutine sf_lresp_ESA + + SUBROUTINE optrot(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci,no,nv, + . xm,ym,zm,xmass) + use commonresp + use commonlogicals + use omp_lib + IMPLICIT NONE + + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci + integer ::io1,io2,iv1,iv2,iwrk,jwrk + integer ::maxconf,moci,no,nv + integer ::iconf(maxconf,2) + integer*8 ::lin8 + + real*8 ::xl(moci*(moci+1)/2) + real*8 ::yl(moci*(moci+1)/2) + real*8 ::zl(moci*(moci+1)/2) + + real*8 ::xm(moci*(moci+1)/2) + real*8 ::ym(moci*(moci+1)/2) + real*8 ::zm(moci*(moci+1)/2) + real*8 ::xmass,refindex,vorfaktor,reffaktor,refval + logical::da + + real*4 ::apb(nci*(nci+1)/2) + real*4 ::amb(nci*(nci+1)/2) + real*4, allocatable ::inv_amb(:) + real*4, allocatable ::inv_resp(:) + real*8, allocatable ::XpY(:,:),XmY(:,:) + real*4 ::omega + real*4, allocatable ::freq(:) + real*8 ::alpha_xx,alpha_xy,alpha_xz + real*8 ::alpha_yy,alpha_yz + real*8 ::alpha_zz + character*1 ::uplo + integer ::info + integer, allocatable ::ipiv(:) + real*4, allocatable ::work (:) + + integer ::ix,iy,iz + + real*4 ::start_time,end_time + + character(len=14):: dummy + + logical :: file_exists + + integer :: nlines + + write(*,*) + write(*,*)'==================================================== + .==================' + write(*,*)' Welcome in linear response sTD-DFT + . program' + write(*,*)'==================================================== + .==================' + + INQUIRE(FILE="wavelength", EXIST=file_exists) + if(file_exists==.false.)then + num_freq=12 + allocate(freq(12)) + freq(1)=45.56335/1900.00 + freq(2)=45.56335/1500.00 + freq(3)=45.56335/1064.00 + freq(4)=45.56335/929.00 + freq(5)=45.56335/794.00 + freq(6)=45.56335/713.00 + freq(7)=45.56335/632.80 + freq(8)=45.56335/589.30 + freq(9)=45.56335/579.00 + freq(10)=45.56335/546.00 + freq(11)=45.56335/436.00 + freq(12)=45.56335/365.00 + else + nlines = 0 + open(unit=101,file='wavelength',form='formatted') + DO + READ (101,*, END=10) + nlines = nlines + 1 + END DO + 10 rewind (101) + num_freq=nlines + allocate(freq(num_freq)) + Do i=1, num_freq + read(101,*)freq(i) + write(*,*)freq(i) + freq(i)=45.56335/freq(i) + enddo + close(101) + endif + + write(*,*) 'Optical rotation alpha[grad*cm^3*g^-1*dm^-1]' + write(*,*) 'including Lorentz factor for common solvent (n=1.4)' + write(*,*) 'lambda [eV] alpha alpha(n=1.0) Phi + . Phi(n=1.0)' + ! shift + + +c refractive index of solvent + refindex=1.4d0 + reffaktor=(refindex**2+2.0d0)/3.0d0 + + allocate(inv_amb(nci*(nci+1)/2)) + inv_amb=amb + uplo='U' + allocate(ipiv(1:nci),work(1:nci)) + call ssptrf(uplo,nci,inv_amb,ipiv,info) + call ssptri(uplo,nci,inv_amb,ipiv,work,info) + deallocate(ipiv,work) + + allocate( XpY(nci,3)) + allocate( XmY(nci,3)) + + Do ii=1, num_freq + omega=freq(ii) + XpY(:,:)=0.0 + call cpu_time(start_time) + allocate(inv_resp(nci*(nci+1)/2)) + inv_resp=apb-omega**2.0*inv_amb + + + uplo='U' + allocate(ipiv(1:nci),work(1:nci)) + call ssptrf(uplo,nci,inv_resp,ipiv,info) + call ssptri(uplo,nci,inv_resp,ipiv,work,info) + deallocate(ipiv,work) + + + +!$omp parallel private(i,j,jk,io,iv,idum1,idum2,ij) +!$omp& reduction (+:XpY) +!$omp do + Do i=1, nci + Do j=1, nci + jk=lin8(i,j) + io=iconf(j,1) + iv=iconf(j,2) + idum1=max(io,iv) + idum2=min(io,iv) + ij=idum2+idum1*(idum1-1)/2 + XpY(i,1)=XpY(i,1)-2.0*xl(ij)*dble(inv_resp(jk)) + XpY(i,2)=XpY(i,2)-2.0*yl(ij)*dble(inv_resp(jk)) + XpY(i,3)=XpY(i,3)-2.0*zl(ij)*dble(inv_resp(jk)) + enddo + enddo +!$omp end do +!$omp end parallel + + !(X-Y)=omega*(A-B)^-1 (X+Y) + XmY=0.0 +!$omp parallel private(i,j,ij,ix) +!$omp& reduction (+:XmY) +!$omp do + Do i=1,nci + Do j=1,nci + ij=lin8(i,j) + Do ix=1,3 + XmY(i,ix)=XmY(i,ix)+ + . dble(inv_amb(ij))*XpY(j,ix)!*dble(omega) !because beta=-Tr(G)/omega/3 + enddo + enddo + enddo +!$omp end do +!$omp end parallel + if(nto)then + write(dummy,'(a,i0)')'1-',ii + open(unit=14,file=dummy) + write(14,*)XmY(:,1)*dble(omega) + close(14) + write(dummy,'(a,i0)')'2-',ii + open(unit=14,file=dummy) + write(14,*)XmY(:,2)*dble(omega) + close(14) + write(dummy,'(a,i0)')'3-',ii + open(unit=14,file=dummy) + write(14,*)XmY(:,3)*dble(omega) + close(14) + endif + alpha_xx=0.0 + alpha_yy=0.0 + alpha_zz=0.0 +!$omp parallel private(i,io,iv,idum1,idum2,ij) +!$omp& reduction(+:alpha_xx,alpha_yy,alpha_zz) +!$omp do + Do j=1, nci + io=iconf(j,1) + iv=iconf(j,2) + idum1=max(io,iv) + idum2=min(io,iv) + ij=idum2+idum1*(idum1-1)/2 + alpha_xx=alpha_xx+(xm(ij)*XmY(j,1))*2.0 ! minus included + alpha_yy=alpha_yy+(ym(ij)*XmY(j,2))*2.0 + alpha_zz=alpha_zz+(zm(ij)*XmY(j,3))*2.0 + enddo +!$omp end do +!$omp end parallel + alpha_xx=alpha_xx/3.0 + alpha_yy=alpha_yy/3.0 + alpha_zz=alpha_zz/3.0 + !G_xx= omega*-alpha_xx + !beta=-Tr(G)/omega/3=Tr(alpha)/3 + + vorfaktor=(38652./xmass)*(589.3/(45.56335/omega))**2 + if(ii==8)then + refval=0 + inquire(file='.ref',exist=da) + if(da)then + open(unit=33,file='.ref') + read(33,*)refval + close(33) + endif + write(*,142)45.56335/omega,omega*27.21139, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*reffaktor*vorfaktor, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*vorfaktor, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*reffaktor*vorfaktor*xmass/100.0, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*vorfaktor*xmass/100.0, + .refval,' Na D-line ' + else + write(*,143)45.56335/omega,omega*27.21139, + .237.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*reffaktor*vorfaktor, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*vorfaktor, + .237.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*reffaktor*vorfaktor*xmass/100.0, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*vorfaktor*xmass/100.0 + endif + + deallocate(inv_resp) + enddo + deallocate(XpY,XmY) + deallocate(inv_amb) +111 format(A15,F20.6) +1111 format(A22,2A20) +2222 format(A2,3F20.6) +3333 format(A16,F7.1,A1,F7.1,A1) + 142 format(f6.1,f6.2,5f12.2,a11) + 143 format(f6.1,f6.2,4f12.2) + end subroutine optrot + + SUBROUTINE optrot_velo(nci,apb,amb,iconf,maxconf,xv,yv,zv,moci, + . no,nv, + . xm,ym,zm,xmass) + use commonresp + use commonlogicals + use omp_lib + IMPLICIT NONE + + integer ::i,j,k,ii,jj,kk,ij,jk,ab,io,iv,idum1,idum2,nci + integer ::io1,io2,iv1,iv2,iwrk,jwrk + integer ::maxconf,moci,no,nv + integer ::iconf(maxconf,2) + integer*8 ::lin8 + + real*8 ::xv(moci*(moci+1)/2) + real*8 ::yv(moci*(moci+1)/2) + real*8 ::zv(moci*(moci+1)/2) + + real*8 ::xvelo(moci*(moci+1)/2) + real*8 ::yvelo(moci*(moci+1)/2) + real*8 ::zvelo(moci*(moci+1)/2) + + real*8 ::xm(moci*(moci+1)/2) + real*8 ::ym(moci*(moci+1)/2) + real*8 ::zm(moci*(moci+1)/2) + real*8 ::xmass,refindex,vorfaktor,reffaktor,refval + logical::da + + real*4 ::apb(nci*(nci+1)/2) + real*4 ::amb(nci*(nci+1)/2) + real*4, allocatable ::inv_amb(:) + real*4, allocatable ::inv_resp(:) + real*8, allocatable ::XpY(:,:),XmY(:,:) + real*4 ::omega + real*4, allocatable ::freq(:) + real*8 ::alpha_xx,alpha_xy,alpha_xz + real*8 ::alpha_yy,alpha_yz + real*8 ::alpha_zz + character*1 ::uplo + integer ::info + integer, allocatable ::ipiv(:) + real*4, allocatable ::work (:) + + integer ::ix,iy,iz + + real*4 ::start_time,end_time + + character(len=14):: dummy + + logical :: file_exists + + integer :: nlines + + write(*,*) + write(*,*)'==================================================== + .==================' + write(*,*)' Welcome in linear response sTD-DFT + . program' + write(*,*)'==================================================== + .==================' + write(*,*)'VELOCITY REPRESENTATION (ORIGIN INDEPENDENT)' + INQUIRE(FILE="wavelength", EXIST=file_exists) + if(file_exists==.false.)then + num_freq=12 + allocate(freq(12)) + freq(1)=45.56335/1900.00 + freq(2)=45.56335/1500.00 + freq(3)=45.56335/1064.00 + freq(4)=45.56335/929.00 + freq(5)=45.56335/794.00 + freq(6)=45.56335/713.00 + freq(7)=45.56335/632.80 + freq(8)=45.56335/589.30 + freq(9)=45.56335/579.00 + freq(10)=45.56335/546.00 + freq(11)=45.56335/436.00 + freq(12)=45.56335/365.00 + else + nlines = 0 + open(unit=101,file='wavelength',form='formatted') + DO + READ (101,*, END=10) + nlines = nlines + 1 + END DO + 10 rewind (101) + num_freq=nlines + allocate(freq(num_freq)) + Do i=1, num_freq + read(101,*)freq(i) + write(*,*)freq(i) + freq(i)=45.56335/freq(i) + enddo + close(101) + endif + write(*,*) 'Optical rotation alpha[grad*cm^3*g^-1*dm^-1]' + write(*,*) 'including Lorentz factor for common solvent (n=1.4)' + write(*,*) 'lambda [eV] alpha alpha(n=1.0) Phi + . Phi(n=1.0)' + +c refractive index of solvent + refindex=1.4d0 + reffaktor=(refindex**2+2.0d0)/3.0d0 + + allocate(inv_amb(nci*(nci+1)/2)) + inv_amb=amb + uplo='U' + allocate(ipiv(1:nci),work(1:nci)) + call ssptrf(uplo,nci,inv_amb,ipiv,info) + call ssptri(uplo,nci,inv_amb,ipiv,work,info) + deallocate(ipiv,work) + + ! mu_v= -(A-B)**(-1)*nabla + + xvelo=0.0 + yvelo=0.0 + zvelo=0.0 +!$omp parallel private(i,j,jk,io,iv,idum1,idum2,ij,ab) +!$omp& reduction (+:xvelo,yvelo,zvelo) +!$omp do + Do i=1, nci + io=iconf(i,1) + iv=iconf(i,2) + idum1=max(io,iv) + idum2=min(io,iv) + ab=idum2+idum1*(idum1-1)/2 + Do j=1, nci + jk=lin8(i,j) + io=iconf(j,1) + iv=iconf(j,2) + idum1=max(io,iv) + idum2=min(io,iv) + ij=idum2+idum1*(idum1-1)/2 + xvelo(ab)=xvelo(ab)+xv(ij)*dble(inv_amb(jk)) + yvelo(ab)=yvelo(ab)+yv(ij)*dble(inv_amb(jk)) + zvelo(ab)=zvelo(ab)+zv(ij)*dble(inv_amb(jk)) + enddo + enddo +!$omp end do +!$omp end parallel + + allocate( XpY(nci,3)) + allocate( XmY(nci,3)) + + Do ii=1, num_freq + omega=freq(ii) + XpY(:,:)=0.0 + call cpu_time(start_time) + allocate(inv_resp(nci*(nci+1)/2)) + inv_resp=apb-omega**2.0*inv_amb + + + uplo='U' + allocate(ipiv(1:nci),work(1:nci)) + call ssptrf(uplo,nci,inv_resp,ipiv,info) + call ssptri(uplo,nci,inv_resp,ipiv,work,info) + deallocate(ipiv,work) + + + +!$omp parallel private(i,j,jk,io,iv,idum1,idum2,ij) +!$omp& reduction (+:XpY) +!$omp do + Do i=1, nci + Do j=1, nci + jk=lin8(i,j) + io=iconf(j,1) + iv=iconf(j,2) + idum1=max(io,iv) + idum2=min(io,iv) + ij=idum2+idum1*(idum1-1)/2 + XpY(i,1)=XpY(i,1)-2.0*xvelo(ij)*dble(inv_resp(jk)) + XpY(i,2)=XpY(i,2)-2.0*yvelo(ij)*dble(inv_resp(jk)) + XpY(i,3)=XpY(i,3)-2.0*zvelo(ij)*dble(inv_resp(jk)) + enddo + enddo +!$omp end do +!$omp end parallel + + !(X-Y)=omega*(A-B)^-1 (X+Y) + XmY=0.0 +!$omp parallel private(i,j,ij,ix) +!$omp& reduction (+:XmY) +!$omp do + Do i=1,nci + Do j=1,nci + ij=lin8(i,j) + Do ix=1,3 + XmY(i,ix)=XmY(i,ix)+ + . dble(inv_amb(ij))*XpY(j,ix)!*dble(omega) because beta=-Tr(G)/omega/3 + enddo + enddo + enddo +!$omp end do +!$omp end parallel + if(nto)then + write(dummy,'(a,i0)')'1-',ii + open(unit=14,file=dummy) + write(14,*)XmY(:,1)*dble(omega) + close(14) + write(dummy,'(a,i0)')'2-',ii + open(unit=14,file=dummy) + write(14,*)XmY(:,2)*dble(omega) + close(14) + write(dummy,'(a,i0)')'3-',ii + open(unit=14,file=dummy) + write(14,*)XmY(:,3)*dble(omega) + close(14) + endif + alpha_xx=0.0 + alpha_yy=0.0 + alpha_zz=0.0 +!$omp parallel private(i,io,iv,idum1,idum2,ij) +!$omp& reduction(+:alpha_xx,alpha_yy,alpha_zz) +!$omp do + Do j=1, nci + io=iconf(j,1) + iv=iconf(j,2) + idum1=max(io,iv) + idum2=min(io,iv) + ij=idum2+idum1*(idum1-1)/2 + alpha_xx=alpha_xx+(xm(ij)*XmY(j,1))*2.0 ! minus included + alpha_yy=alpha_yy+(ym(ij)*XmY(j,2))*2.0 + alpha_zz=alpha_zz+(zm(ij)*XmY(j,3))*2.0 + enddo +!$omp end do +!$omp end parallel + alpha_xx=alpha_xx/3.0 + alpha_yy=alpha_yy/3.0 + alpha_zz=alpha_zz/3.0 + !G_xx= omega*-alpha_xx + !beta=-Tr(G)/omega/3=Tr(alpha)/3 + + vorfaktor=(38652./xmass)*(589.3/(45.56335/omega))**2 + if(ii==8)then + refval=0 + inquire(file='.ref',exist=da) + if(da)then + open(unit=33,file='.ref') + read(33,*)refval + close(33) + endif + write(*,142)45.56335/omega,omega*27.21139, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*reffaktor*vorfaktor, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*vorfaktor, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*reffaktor*vorfaktor*xmass/100.0, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*vorfaktor*xmass/100.0, + .refval,' Na D-line ' + else + write(*,143)45.56335/omega,omega*27.21139, + .237.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*reffaktor*vorfaktor, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*vorfaktor, + .237.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*reffaktor*vorfaktor*xmass/100.0, + .235.722*(alpha_xx+alpha_yy+alpha_zz) + ./64604.8*(2.0*137.036/3.0)*vorfaktor*xmass/100.0 + endif + + deallocate(inv_resp) + enddo + deallocate(XpY,XmY) + deallocate(inv_amb) +111 format(A15,F20.6) +1111 format(A22,2A20) +2222 format(A2,3F20.6) +3333 format(A16,F7.1,A1,F7.1,A1) + 142 format(f6.1,f6.2,5f12.2,a11) + 143 format(f6.1,f6.2,4f12.2) + end subroutine optrot_velo diff --git a/main.f b/main.f index cd5701f..de6df7e 100644 --- a/main.f +++ b/main.f @@ -14,8 +14,8 @@ ! ! You should have received a copy of the GNU Lesser General Public License ! along with stda. If not, see . -! - program acis_prog +! + program acis_prog use stdacommon ! mostly input and primitive data use kshiftcommon ! kshiftvariables use commonlogicals @@ -27,9 +27,9 @@ program acis_prog integer, allocatable :: ccspin(:) real*8, allocatable ::xyz(:,:) real*8 xx(10),alpha,beta,ptlim - character*79 dummy + character*79 dummy character*79 fname - character*8 method + character*8 method integer imethod,inpchk,mform,nvec logical molden,da,chkinp,xtbinp @@ -55,7 +55,7 @@ program acis_prog write(*,'(a,a)')'===============================================', . '=======================' write(*,*) -c defaults +c defaults ptlim=1.7976931348623157d308 ! energy range that will be scanned by PT (we use just a large number) thre=7.0 ! energy range for primary CSF space @@ -63,13 +63,13 @@ program acis_prog beta=-100.0d0 ! otherwise global hybrid defaults will be used c the following value yields converged UV spectra for several members of -c of the DYE12 set +c of the DYE12 set thrp=1.d-4 mform=1 ! mform is the "style" specifier for Molden input, by default try TM input: ORCA/XTB = 0, TM=1,Molpro=2, Terachem/Gaussian=3 - rpachk=.false. ! sTD-DFT wanted? - triplet=.false. ! triplet excitations wanted? + rpachk=.false. ! sTD-DFT wanted? + triplet=.false. ! triplet excitations wanted? eigvec=.false. ! eigenvector printout wanted? nvec=0 ! if so, how many vecs? printexciton=.false. ! print information for exciton coupling program @@ -80,15 +80,15 @@ program acis_prog fname='molden.input' xtbinp=.false. !use xtbinput? screen=.false. ! prescreen configurations (by Schwarz-type screening) - + ! Kia shifting defaults dokshift=.false. shftmax=0.50d0 ! maximum Kia shift in eV shftwidth=0.10d0 ! damping threshold in eV shftsteep=4.0d0 ! steepness -c read the tm2xx file, otherwise (-f option) the tm2molden file - molden=.true. +c read the tm2xx file, otherwise (-f option) the tm2molden file + molden=.true. ax=-1 imethod=0 inpchk=0 @@ -99,16 +99,18 @@ program acis_prog smp2=.false. bfw=.false. spinflip=.false. + sf_s2=.false. rw=.false. pt_off=.false. + optrot=.false. ! check for input file inquire(file='.STDA',exist=da) if(da)then call readinp(ax,thre,alpha,beta) - endif + endif - do i=1,command_argument_count() + do i=1,command_argument_count() call getarg(i,dummy) if(index(dummy,'-fo').ne.0)then call getarg(i+1,fname) @@ -120,7 +122,7 @@ program acis_prog molden=.true. inpchk=inpchk+1 endif - + if(index(dummy,'-ax').ne.0)then ! get amout of Fock-exchange call getarg(i+1,dummy) call readl(79,dummy,xx,nn) @@ -158,7 +160,7 @@ program acis_prog call readl(79,dummy,xx,nn) if(nn.gt.0) ptlim=xx(1) endif - + if(index(dummy,'-xtb').ne.0) then ! xTB input, set defaults ! two other dirty paramters (not ! really fitted in function @@ -192,7 +194,7 @@ program acis_prog call readl(79,dummy,xx,nn) num_freq=int(xx(1)) endif - + if(index(dummy,'-aresp').ne.0) then aresp=.true. rpachk=.true. @@ -200,7 +202,7 @@ program acis_prog call readl(79,dummy,xx,nn) num_freq=int(xx(1)) endif - + if(index(dummy,'-2PA').ne.0)then ! Do response function TPA=.true. rpachk=.true. @@ -208,7 +210,7 @@ program acis_prog call readl(79,dummy,xx,nn) num_trans=int(xx(1)) endif - + if(index(dummy,'-s2s').ne.0)then ! Do response function ESA=.true. !rpachk=.true. @@ -216,22 +218,22 @@ program acis_prog call readl(79,dummy,xx,nn) state2opt=int(xx(1)) endif - + ! if(index(dummy,'-MP2').ne.0)then ! Do mp2 ! smp2=.true. ! endif - + if(index(dummy,'-rw').ne.0)then ! Do mp2 rw=.true. call getarg(i+1,dummy) call readl(79,dummy,xx,nn) if(xx(1)==1) pt_off=.true. endif - + if(index(dummy,'-BFW').ne.0)then bfw=.true. endif - + if(index(dummy,'-sf').ne.0)then ! Do spinflip spinflip=.true. if(xtbinp) then @@ -239,26 +241,38 @@ program acis_prog ax=0.36d0 endif endif - + if(index(dummy,'-spin').ne.0)then + sf_s2=.true. + endif + + if(index(dummy,'-oprot').ne.0)then + rpachk=.true. + optrota=.true. + velo_OR=.false. + call getarg(i+1,dummy) + call readl(79,dummy,xx,nn) + if(xx(1)==1)velo_OR=.true. + endif + if(index(dummy,'-nto').ne.0)then ! Do nto nto=.true. call getarg(i+1,dummy) call readl(79,dummy,xx,nn) Nnto=int(xx(1)) endif - + if(index(dummy,'-chk').ne.0) chkinp=.true. ! do input check - if(index(dummy,'-vectm').ne.0)then + if(index(dummy,'-vectm').ne.0)then eigvec=.true. ! print eigenvectors call getarg(i+1,dummy) call readl(79,dummy,xx,nn) if(nn.gt.0) nvec=dnint(xx(1)) endif ! print transition dipole moments - if(index(dummy,'-excprint').ne.0) printexciton=.true. - if(index(dummy,'-oldtda').ne.0) velcorr=.false. + if(index(dummy,'-excprint').ne.0) printexciton=.true. + if(index(dummy,'-oldtda').ne.0) velcorr=.false. if(index(dummy,'-aniso').ne.0) aniso=.true. - enddo + enddo ! if(ptlim.gt.1.0d308) ptlim = 3.0d0*thre ! set ptlimit to a multiple of thre @@ -266,18 +280,18 @@ program acis_prog ccccccccccccccccccccccccccccccc c check the input -ccccccccccccccccccccccccccccccc +ccccccccccccccccccccccccccccccc if(inpchk.gt.1) stop 'please specify only one input file!' - if(inpchk.lt.1) stop 'no input file specified!' + if(inpchk.lt.1) stop 'no input file specified!' if(molden) then write(*,*) 'reading a molden input...' - else if (xtbinp) then + else if (xtbinp) then write(*,*) 'reading an xTB output...' - else + else write(*,*) 'reading a tm2stda file...' chkinp=.false. endif - if(ax.lt.0.and..not.chkinp) then + if(ax.lt.0.and..not.chkinp) then stop 'specify Fock exchange via -ax ' endif @@ -297,15 +311,15 @@ program acis_prog c compare input and inpchk if(inpchk.ne.imethod)stop'U/R-KS input doesnot match moldenfile' - else if(xtbinp) then ! read parameters from xTB input + else if(xtbinp) then ! read parameters from xTB input call readxtb0(imethod,ncent,nmo,nbf,nprims) else - call readbas0a(0,ncent,nmo,nbf,nprims,fname) + call readbas0a(0,ncent,nmo,nbf,nprims,fname) endif if(nprims.eq.0.or.ncent.eq.0.or.nmo.eq.0) - .stop 'read error' + .stop 'read error' ccccccccccccccccccccccccccccccc c allocate mo vectors @@ -349,8 +363,8 @@ program acis_prog call readmold(mform,imethod,ncent,nmo,nbf,nprims,cc, . ccspin,icdim,fname) else if(xtbinp) then - call readxtb(imethod,ncent,nmo,nbf,nprims,cc) - if(imethod.eq.2) then + call readxtb(imethod,ncent,nmo,nbf,nprims,cc) + if(imethod.eq.2) then do i=1,nmo/2 ccspin(i)=1 enddo @@ -360,21 +374,21 @@ program acis_prog endif else if(imethod.eq.1) call readbasa(1,imethod,ncent,nmo,nbf,nprims,cc, - .icdim,fname,iaobas) + .icdim,fname,iaobas) if(imethod.eq.2) call readbasb(1,imethod,ncent,nmo,nbf,nprims,cc, - .ccspin,icdim,fname,iaobas) + .ccspin,icdim,fname,iaobas) endif if(imethod.eq.1) deallocate( ccspin ) ccccccccccccccccccccccccccccccccc c precalculate primitive data ccccccccccccccccccccccccccccccccc - call intslvm(ncent,nmo,nbf,nprims) + call intslvm(ncent,nmo,nbf,nprims) nao=nbf if(nao.eq.0)nao=nprims - - - + + + if(nto)then ! ! Have everything necessary to compute nice NTOs @@ -411,12 +425,12 @@ program acis_prog open(unit=12,file='fnorm') close(12,status='delete') endif - - 21 format(3i10,3x,2f16.9) - - 11 format(3i5,2f9.4) + + 21 format(3i10,3x,2f16.9) + + 11 format(3i5,2f9.4) deallocate(ipat,ipty,ipao,exip,cxip,atnam,eta) - + cccccccccccccccccccccccccccccccccccccc c if input check (-chk) is used c @@ -430,22 +444,22 @@ program acis_prog close(36,status='delete') deallocate( cc ) if(imethod.eq.2) deallocate( ccspin ) - if(mform.gt.3) stop 'unreadable format' + if(mform.gt.3) stop 'unreadable format' write(*,*)'Restarting and trying different input style...' deallocate(iaoat,occ,eps,co) goto 888 else call inputcheck_printout(mform,.true.,imethod,nbf,nmo) call exit(0) ! leave program - endif + endif else - if(.not.xtbinp) then + if(.not.xtbinp) then call inputcheck_printout(mform,.false.,imethod,nbf,nmo) endif endif cccccccccccccccccccccccccccccccccccccccccc c make an additional printout whether c -c restricted or unrestricted calculation c +c restricted or unrestricted calculation c c will be done c cccccccccccccccccccccccccccccccccccccccccc write(*,*)' ' @@ -466,7 +480,7 @@ program acis_prog endif write(*,*) ' ' ccccccccccccccccccccccccccccccccc -c stda +c stda ccccccccccccccccccccccccccccccccc allocate(xyz(4,ncent),stat=ierr) if(ierr.ne.0) stop 'allocation failed in main for xyz' @@ -475,17 +489,17 @@ program acis_prog xyz(1:4,i)=co(i,1:4) enddo deallocate(co) - + ! if(smp2)then ! call sMP(ncent,nmo,nao,xyz,cc,eps,occ,iaoat, ! . ax,alpha,beta) ! endif - + if(spinflip)then call sfstda(ncent,nmo,nao,xyz,cc,eps,occ,ccspin,iaoat,thre, . thrp,ax,alpha,beta,ptlim,nvec) endif - + if (imethod.eq.1) then if(rw)then call stda_rw(ncent,nmo,nao,xyz,cc,eps,occ,iaoat,thre, @@ -496,7 +510,7 @@ program acis_prog endif else call sutda(ncent,nmo,nao,xyz,cc,eps,occ,ccspin,iaoat,thre, - . thrp,ax,alpha,beta,ptlim,nvec) + . thrp,ax,alpha,beta,ptlim,nvec) endif call system('date') @@ -526,7 +540,7 @@ subroutine byteout(s,mem) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! this subroutine performs a mulliken population analysis to check if the +! this subroutine performs a mulliken population analysis to check if the ! input data is correct subroutine mulpopcheck(nbf,nmo,c,occ,ireturn) implicit none @@ -537,13 +551,13 @@ subroutine mulpopcheck(nbf,nmo,c,occ,ireturn) integer, intent( out ) :: ireturn real*8, allocatable ::pmul(:,:),ovvec(:),ovmat(:,:),dmat(:,:) allocate( dmat(nbf,nbf), ovmat(nbf,nbf), ovvec(nbf*(nbf+1)/2), - . pmul(nbf,nbf) ) + . pmul(nbf,nbf) ) ireturn=0 write(*,*) ' ' call header('I N P U T C H E C K',0) ! for debugging ! k=0 -! do i=1,nmo +! do i=1,nmo ! write(*,*) 'MO', i ! do j=1,nbf ! k=k+1 @@ -568,8 +582,8 @@ subroutine mulpopcheck(nbf,nmo,c,occ,ireturn) dmat(i,j)=dmat(i,j)+occ(k)*c(nbf*(k-1)+i)*c(nbf*(k-1)+j) dmat(j,i)=dmat(i,j) enddo - enddo - enddo + enddo + enddo ! call prmat(6,dmat,nbf,nbf,'density matrix') summe2=0.0d0 ! get number of electrons in the system (reference) @@ -584,13 +598,13 @@ subroutine mulpopcheck(nbf,nmo,c,occ,ireturn) enddo ! call prmat(6,pmul,nbf,nbf,'P * S') deallocate(dmat,ovmat,pmul) ! free memory - if(dabs(summe1-summe2).gt.nmo*1.d-5)then ! since Gaussian comes with accuracy up to 5 decimal points, take a system-dependent accuracy + if(dabs(summe1-summe2).gt.nmo*1.d-5)then ! since Gaussian comes with accuracy up to 5 decimal points, take a system-dependent accuracy write(*,*) 'Number of electrons from Mulliken population' write(*,'(a,x,f18.6,x,f18.6)') 'does not match:',summe1,summe2 write(*,*) 'WARNING: Input style is incorrect!' write(*,*) ' ' ireturn=1 - return + return endif ! print out which input format was detected if(ireturn.eq.0)then @@ -610,7 +624,7 @@ subroutine inputcheck_printout (ires,typ,imethod,nbf,nmo) if (typ) then write(*,*) ' ' write(*,*) ' --- S U C C E S S --- ' - write(*,*) ' ' + write(*,*) ' ' select case(ires) case(0) write(*,*) ' GTO/MO data matches ORCA style' @@ -620,14 +634,14 @@ subroutine inputcheck_printout (ires,typ,imethod,nbf,nmo) write(*,*) ' GTO/MO data matches TURBOMOLE style' write(*,*) ' use "-sty 1" flag in the future' else - write(*,*) ' GTO/MO data matches Cartesian basis style' + write(*,*) ' GTO/MO data matches Cartesian basis style' write(*,*) ' (compatible with TURBOMOLE/MOLPRO/G09)' write(*,*) ' use the following flag in the future' write(*,*) ' TURBOMOLE: "-sty 1" ' write(*,*) ' GAUSSIAN 09: "-sty 1" ' write(*,*) ' MOLPRO: "-sty 2" ' endif - case(2) + case(2) write(*,*) ' GTO/MO data matches MOLPRO style' write(*,*) ' use "-sty 2" flag in the future' case(3) @@ -656,8 +670,8 @@ subroutine inputcheck_printout (ires,typ,imethod,nbf,nmo) open(unit=39,file='zvint',form='unformatted',status='old') close(39,status='delete') return - else - ! standard run: no check + else + ! standard run: no check write(*,*) ' ' write(*,*) 'Skipping input check...' select case(ires) @@ -669,7 +683,7 @@ subroutine inputcheck_printout (ires,typ,imethod,nbf,nmo) else write(*,*) 'Assuming TM/MOLPRO/G09 style (-sty 1)' endif - case(2) + case(2) write(*,*) 'Assuming MOLPRO style (-sty 2)' case(3) write(*,*) 'Assuming TERACHEM style (-sty 3)' diff --git a/sfstda.f b/sfstda.f index 0e3e326..943cf03 100644 --- a/sfstda.f +++ b/sfstda.f @@ -19,10 +19,10 @@ ! written by Marc de Wegifosse 2018-2019 SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, - . ax,alphak,betaj,fthr,nvec) + . ax,alphak,betaj,fthr,nvec) use commonlogicals use omp_lib - IMPLICIT NONE + IMPLICIT NONE c input: integer ncent,nmo,nao integer iaoat(*) @@ -34,8 +34,8 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, c local variables: integer :: STATUS - -c l=dipole lengths, v=dipole velocity (d/dxyz), m=angular momentum + +c l=dipole lengths, v=dipole velocity (d/dxyz), m=angular momentum real*8, allocatable ::xla(:),yla(:),zla(:) ! real*8, allocatable ::xva(:),yva(:),zva(:) ! real*8, allocatable ::xma(:),yma(:),zma(:) @@ -43,15 +43,15 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ! real*8, allocatable ::xvb(:),yvb(:),zvb(:) ! real*8, allocatable ::xmb(:),ymb(:),zmb(:) real*8, allocatable ::help(:),scr(:),dum(:),x(:,:) -c MOs, orbital energies and CSF printout stuff +c MOs, orbital energies and CSF printout stuff real*8, allocatable ::cca(:),epsia(:) real*8, allocatable ::ccb(:),epsib(:) real*8, allocatable ::umerk(:,:) real*8, allocatable ::umrkx(:,:),umrky(:,:),umrkz(:,:) real*8, allocatable ::rvp(:) -c stuff for diag of TDA matrix -c critical regarding memory +c stuff for diag of TDA matrix +c critical regarding memory integer info,lwork,liwork,il,iu,nfound real*4, allocatable ::uci (:,:) real*4, allocatable ::eci (:) @@ -62,7 +62,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, integer,allocatable::isuppz(:) c Löwdin MOs, repulsion terms, charges and half-transformed stuff -c critical regarding memory +c critical regarding memory real*8, allocatable ::clowa(:) real*8, allocatable ::clowb(:) real*4, allocatable ::gamj(:,:) @@ -76,7 +76,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, real*4 sdot c the maximum size of the TDA expansion space -! integer maxconf1,maxconf2 +! integer maxconf1,maxconf2 ! parameter (maxconf1=500000) !c the maximum size of the TDA pt2 space ! parameter (maxconf2=maxconf1*10) @@ -95,7 +95,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, real*4, allocatable ::apb(:) real*4, allocatable ::ambsqr(:) -c intermediates +c intermediates real*8 omax,vmin,pert,de,ek,ej,ak,xc,rabx,ef real*8 pp,hilf,uu,sss,rl,rv,rm,time,coc(3) real*8 fl,fv,ec,p23,xp,umax,xvu,yvu,zvu,xmu,ymu,zmu,xlu,ylu,zlu @@ -117,16 +117,23 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ! variables for vector printout integer, allocatable :: vecchka(:),vecchkb(:) integer nvec,jhomo,jhomoa,jhomob -c atomic Hubbard parameters +c atomic Hubbard parameters real*8 eta(94) -c atomic masses +c atomic masses common /amass/ ams(107) real*8 ams character*79 dummy - -c just a printout +c For the computation of + real*8, allocatable :: i_alpha_j_beta(:,:) + real*8, allocatable :: i_alpha_a_beta(:,:) + real*8 :: sum_ij,S2_UHF + real*8, allocatable :: Xia(:),Xiaib(:),Xiaka(:) + real*8, allocatable :: S2(:) + integer kk + +c just a printout call header('s T D A',0) - call cpu_time(start_time) + call cpu_time(start_time) thr =thr /27.211385050d0 c estimate the orbital energy window which corresponds to the desired c spectra range thr @@ -140,15 +147,15 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, if(occ(i).gt.0.99.and.eps(i).gt.omax) omax=eps(i) if(occ(i).lt.0.01.and.eps(i).lt.vmin) vmin=eps(i) enddo - - + + ! optional: if eigenvectors are wanted in TM format, check now how many occupied there are in general ! this is needed to get the CSF sorting of TM jhomo=0 jhomoa=0 jhomob=0 do i=1,nmo - if(occ(i).gt.0.990d0) then + if(occ(i).gt.0.990d0) then jhomo=jhomo+1 if(csp(i).eq.1) jhomoa = jhomoa + 1 if(csp(i).eq.2) jhomob = jhomob + 1 @@ -191,7 +198,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, . ccb(nao*mocib),epsib(mocib) . ) - write(*,*) 'Active space MOs in TDA : ',moci,mocia,mocib + write(*,*) 'Active space MOs in TDA : ',moci,mocia,mocib mocia=0 mocib=0 @@ -216,13 +223,13 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, endif endif enddo - + no = moci noa = mocia nob = mocib - + do i = 1,nmo - if(occ(i).lt.0.01.and.eps(i).lt.vthr)then + if(occ(i).lt.0.01.and.eps(i).lt.vthr)then moci=moci+1 if(csp(i).eq.1) then mocia = mocia + 1 @@ -245,14 +252,14 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, nva = mocia - noa nvb = mocib - nob - write(*,*) 'Occupied active MOs in TDA: ', no,noa,nob + write(*,*) 'Occupied active MOs in TDA: ', no,noa,nob if(noa<=nob)then write(*,*)'warning number of alpha < beta' write(*,*)'only with alpha > beta' call exit(status) endif - write(*,*) 'Virtual active MOs in TDA: ', nv,nva,nvb - if((noa.eq.0.or.nva.eq.0).and.(nob.eq.0.or.nvb.eq.0)) then + write(*,*) 'Virtual active MOs in TDA: ', nv,nva,nvb + if((noa.eq.0.or.nva.eq.0).and.(nob.eq.0.or.nvb.eq.0)) then stop 'no CSF, increase energy threshold (-e option)' endif @@ -268,7 +275,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, kconfb=0 -c we arrange MOS according to energy from 1:HOMO to LUMO:MOCI +c we arrange MOS according to energy from 1:HOMO to LUMO:MOCI c (in TM they come in irreps) write(*,*) 'Sorting MOs ...' c sort for E diag @@ -310,14 +317,14 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, call onetri(1,help,dum,scr,ccb,nao,mocib) call shrink(mocib,dum,zlb) close(33,status='delete') -! -! +! +! ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! c ! c magnetic dipole ! c ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -! +! open(unit=34,file='xmint',form='unformatted',status='old') ! read(34) help ! call onetri(-1,help,dum,scr,cca,nao,mocia) @@ -325,7 +332,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ! call onetri(-1,help,dum,scr,ccb,nao,mocib) ! call shrink(mocib,dum,xmb) close(34,status='delete') -! +! open(unit=35,file='ymint',form='unformatted',status='old') ! read(35) help ! call onetri(-1,help,dum,scr,cca,nao,mocia) @@ -333,7 +340,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ! call onetri(-1,help,dum,scr,ccb,nao,mocib) ! call shrink(mocib,dum,ymb) close(35,status='delete') -! +! open(unit=36,file='zmint',form='unformatted',status='old') ! read(36) help ! call onetri(-1,help,dum,scr,cca,nao,mocia) @@ -341,13 +348,13 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ! call onetri(-1,help,dum,scr,ccb,nao,mocib) ! call shrink(mocib,dum,zmb) close(36,status='delete') -! +! ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! c ! c velocity dipole ! c ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -! +! open(unit=37,file='xvint',form='unformatted',status='old') ! read(37) help ! call onetri(-1,help,dum,scr,cca,nao,mocia) @@ -355,15 +362,15 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ! call onetri(-1,help,dum,scr,ccb,nao,mocib) ! call shrink(mocib,dum,xvb) close(37,status='delete') -! +! open(unit=38,file='yvint',form='unformatted',status='old') ! read(38) help ! call onetri(-1,help,dum,scr,cca,nao,mocia) ! call shrink(mocia,dum,yva) ! call onetri(-1,help,dum,scr,ccb,nao,mocib) ! call shrink(mocib,dum,yvb) - close(38,status='delete') -! + close(38,status='delete') +! open(unit=39,file='zvint',form='unformatted',status='old') ! read(39) help ! call onetri(-1,help,dum,scr,cca,nao,mocia) @@ -373,18 +380,18 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, close(39,status='delete') CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -c calc S^1/2 and q(GS) +c calc S^1/2 and q(GS) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC open(unit=40,file='sint',form='unformatted',status='old') read(40) help write(*,*) 'ints done.' close(40,status='delete') - + write(*,*) 'S^1/2 ...' call makel(nao,help,x) - + call dgemm('n','n',nao,mocia,nao,1.d0,x,nao,cca,nao,0.d0,scr,nao) do i=1,mocia @@ -414,7 +421,33 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, enddo write(*,*) 'S^1/2 orthogonalized MO coefficients done.' - + if(sf_s2)then + ! for the computation of + allocate(i_alpha_j_beta(noa,nob)) + i_alpha_j_beta=0.0 + sum_ij=0.0 + do i=1,noa + do j=1,nob + do k=1,nao + i_alpha_j_beta(i,j)=i_alpha_j_beta(i,j)+ + . clowa(k+(i-1)*nao)*clowb(k+(j-1)*nao) + enddo + sum_ij=sum_ij+i_alpha_j_beta(i,j)**2.0 + enddo + enddo + S2_UHF=(noa-nob)/2.0*((noa-nob)/2.0+1.0)+nob-sum_ij + write(*,*) 'Unrestricted ground state =',S2_UHF + allocate(i_alpha_a_beta(noa,nvb)) + i_alpha_a_beta=0.0 + do i=1,noa + do j=1,nvb + do k=1,nao + i_alpha_a_beta(i,j)=i_alpha_a_beta(i,j)+ + . clowa(k+(i-1)*nao)*clowb(k+(j+nob-1)*nao) + enddo + enddo + enddo + endif deallocate(scr,dum,help,x) @@ -440,7 +473,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ! write(*,'(10F7.3)') q1b(1:ncent) ! write(*,'(10F7.3)') (q1a(1:ncent)+q1b(1:ncent)) ! write(*,*) -! write(*,'('' # of alpha and beta electrons in TDA:'',2F8.3)') +! write(*,'('' # of alpha and beta electrons in TDA:'',2F8.3)') ! .sum(q1a(1:ncent)),sum(q1b(1:ncent)) ! write(*,*) @@ -454,15 +487,15 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ak = 1.0d0 c the global parameters of the method: - beta1 = 0.3 - beta2 = 1.0 - + beta1 = 0.3 + beta2 = 1.0 + if(betaj.lt.-99.0d0) then ! if no beta parameter was read in betaj=beta1+beta2*ax endif - + write(*,*) write(*,*) 'ax(DF) : ',ax write(*,*) 's^K : ',ak @@ -478,9 +511,9 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, xmolw=0 do i = 1,ncent ii = idint(xyz(4,i)) -c ams is the atomic mass (mol weight for output file) +c ams is the atomic mass (mol weight for output file) xmolw = xmolw+ams(ii) - do j=1,i + do j=1,i jj = idint(xyz(4,j)) xj = 0.50d0*(eta(ii)+eta(jj)) * (ax+ax*0.4) !!!!!!!!!!!!!!!!!!!!!! SF this seems resonable rabx = dsqrt((xyz(1,i)-xyz(1,j))**2 @@ -546,8 +579,8 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, iconfb(k,2)=iv edb(k)=de endif -c for PT - if(de.gt.thr.and.de.lt.fthr)then ! needs to be on if fthr is specified +c for PT + if(de.gt.thr.and.de.lt.fthr)then ! needs to be on if fthr is specified j=j+1 kconfb(j,1)=io kconfb(j,2)=iv @@ -558,7 +591,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, enddo nexb=k nexptb=j - deallocate(uci) + deallocate(uci) nex = nexb nexpt = nexptb @@ -569,12 +602,12 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, c errors and warning if(nex.lt.1) stop 'No CSF, increase energy threshold (-e option)' -! if(nexa.eq.maxconf1.or.nexb.eq.maxconf1) +! if(nexa.eq.maxconf1.or.nexb.eq.maxconf1) ! . stop 'Primary CSF space exceeded. use -e option!' -! if(nexpta.eq.maxconf2.or.nexptb.eq.maxconf2) +! if(nexpta.eq.maxconf2.or.nexptb.eq.maxconf2) ! . write(*,*)'CSF PT2 space exceeded. try -p option!' - -c sort for E diag in each spin manifold + +c sort for E diag in each spin manifold do 142 ii = 2,nexb i = ii - 1 @@ -593,10 +626,10 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, iconfb(i,m)=iconfb(k,m) iconfb(k,m)=ihilf enddo - 142 continue + 142 continue -c just printout - write(*,*) ' ' +c just printout + write(*,*) ' ' write(*,*)'Ordered frontier alpha orbitals:' write(*,*)' eV' j=max(1,noa-10) @@ -609,7 +642,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, write(*,'(i4,F10.3,F8.1)') i,epsia(i)*27.21139 enddo - write(*,*) ' ' + write(*,*) ' ' write(*,*)'Ordered frontier beta orbitals:' write(*,*)' eV' j=max(1,nob-10) @@ -662,7 +695,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, enddo enddo !$omp end do -!$omp end parallel +!$omp end parallel allocate(pija(ncent,ihilf), stat=ierr) if(ierr.ne.0)stop 'allocation failed for (ij| intermediate' pija=0.0 @@ -704,7 +737,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, . maxconfb,iconfb,kconfb, . ak,ax,edb,edptb, . pija,qabb,thrp,newb) - + deallocate(edptb,kconfb) new = newb @@ -712,7 +745,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, c nroot at this point is the number of primary CSF. The c number of roots in the range 0-thr (as desired) is not -c known but will be determined by the diag routine. +c known but will be determined by the diag routine. nroot = nex write(*,*) 'CSF included by PT:',new nexb = nexb + newb @@ -733,7 +766,7 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, ******************************************************************************** c set the TDA matrix up * - write(*,*)'calculating TDA matrix ...' + write(*,*)'calculating TDA matrix ...' ******************************************************************************** call sfstdamat(nci,nexb,ncent,noa,nob,nvb, . maxconfb,iconfb,ax,edb, @@ -750,41 +783,82 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, write(*,'('' estimated time (min) '',f8.2)') . float(nci)**2*float(nroot)/8.d+8/60. -c if LAPACK does not work +c if LAPACK does not work c allocate(eci(nci),uci(nci,nroot), c . stat=ierr) c call sHQRII(hci,nci,nroot,eci,uci) c faster by a factor of 2-3 - lwork =26*nci + lwork =26*nci liwork=10*nci -c we allocate uci with nci (and not with nroot) as safe +c we allocate uci with nci (and not with nroot) as safe c choice (other values gave segfaults) nroot=min(nci,int(1.5*nroot)) allocate(eci(nci),uci(nci,nci),work(lwork), . iwork(liwork),isuppz(nci),stat=ierr) if(ierr.ne.0)stop 'allocation failed for TDA matrix diag' - vl=-10.0 + vl=-10.0 vu=thr call ssyevr('V','V','U',nci,hci,nci,vl,vu,il,iu,1.e-6, . nfound,eci,uci,nci,isuppz, . work,lwork,iwork,liwork,info) nroot=nfound if(nfound.lt.1) stop 'internal error in diag' - + if(sf_s2)then + ! Compute for each SF states + allocate(Xia(nroot),Xiaib(nroot),Xiaka(nroot),S2(nroot)) + Xia=0.0 + Xiaib=0.0 + Xiaka=0.0 + Do i=1,nroot + !X_iAaB + Do j=1,nci + Xia(i)=Xia(i)+dble(uci(j,i))* + . i_alpha_a_beta(iconfb(j,1),iconfb(j,2)-nob) + enddo + Xia(i)=Xia(i)**2.0 + !XjAaB XjAbB + Do j=1,nci + Do kk=1,nci + if(iconfb(kk,1)==iconfb(j,1))then + Do k=1,noa + Xiaib(i)=Xiaib(i)+dble(uci(j,i))*dble(uci(kk,i))* + . i_alpha_a_beta(k,iconfb(j,2)-nob)* + . i_alpha_a_beta(k,iconfb(kk,2)-nob) + enddo + endif + enddo + enddo + !XiAaB XkAbB + Do j=1,nci + Do kk=1,nci + if(iconfb(kk,2)==iconfb(j,2))then + Do k=1,nob + Xiaka(i)=Xiaka(i)+dble(uci(j,i))*dble(uci(kk,i))* + . i_alpha_j_beta(iconfb(j,1),k)* + . i_alpha_j_beta(iconfb(kk,1),k) + enddo + endif + enddo + enddo + S2(i)=nob+1-sum_ij-Xiaib(i)+Xiaka(i)+Xia(i)+ + . 0.25*(noa-nob-2)*(noa-nob) + enddo + deallocate(Xia,Xiaib,Xiaka,i_alpha_j_beta,i_alpha_a_beta) + endif write(*,'(i5,'' SF-roots found, lowest/highest eigenvalue : '', .2F8.3,i4)') nroot,eci(1)*27.21139,eci(nroot)*27.21139,info if(info.gt.0) stop 'diag error (ssyevr)' - + allocate(umerk(14,nroot)) - + ! largest configurations - + Do i=1,nroot - umax=-1 - kmem=1 + umax=-1 + kmem=1 ! first largest do j=1,nci io=iconfb(j,1) @@ -837,32 +911,42 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, umerk(8,i)=jmax umerk(9,i)=dble(uci(jmem,i)) enddo - - 20 format(i5,2f9.5,3(' ',i4,'(a) ->',i4,'(b) ',F6.2)) - + + 20 format(i5,2f9.5,3(' ',i4,'(a) ->',i4,'(b) ',F6.2),f9.5) + 21 format(i5,2f9.5,3(' ',i4,'(a) ->',i4,'(b) ',F6.2)) + write(*,*)'SF states:' + if(sf_s2)then Do i=1, nroot write(*,20)i,eci(i),eci(i)*27.21139, .int(umerk(1:2,i)),umerk(3,i), .int(umerk(4:5,i)),umerk(6,i), + .int(umerk(7:8,i)),umerk(9,i),S2(i) + enddo + else + Do i=1, nroot + write(*,21)i,eci(i),eci(i)*27.21139, + .int(umerk(1:2,i)),umerk(3,i), + .int(umerk(4:5,i)),umerk(6,i), .int(umerk(7:8,i)),umerk(9,i) enddo - + endif + call sf_lresp_ESA(nci,iconfb,maxconfb,xla,yla,zla,mocia, . xlb,ylb,zlb,mocib, - . noa,nva,nob,nvb,eci,uci,nroot,xmolw,thr) - - + . noa,nva,nob,nvb,eci,uci,nroot,xmolw,thr) + + if(nto)then - + call SFprint_nto(uci,cca,mocia,nci,nroot,nao, . noa,nva,ccb,mocib,iconfb, . maxconfb,nob,nvb) - endif - - + endif + + deallocate(hci) @@ -873,18 +957,18 @@ SUBROUTINE sfstda(ncent,nmo,nao,xyz,c,eps,occ,csp,iaoat,thr,thrp, 13 format( 41x,f8.3,13x,f8.3) 14 format(3x,F6.2,'(',i4,'a','->',i4,'a',')') 15 format(3x,F6.2,'(',i4,'b','->',i4,'b',')') - 16 format(i5,f6.2,f8.1, 5x,i4,' ->',i4,5x,'gap,J:',2f8.3, - . 3x,'Kshft:',f8.3,2x,'locality:',f6.3) + 16 format(i5,f6.2,f8.1, 5x,i4,' ->',i4,5x,'gap,J:',2f8.3, + . 3x,'Kshft:',f8.3,2x,'locality:',f6.3) call cpu_time(end_time) print '("sTD-DFT Time = ",f12.2," minutes.")' . ,(end_time-stda_time)/60.0 write(*,*) 'SF-sTDA done.' - call EXIT(STATUS) + call EXIT(STATUS) return end SUBROUTINE sfstda - - + + subroutine sfstdamat(nci,nexb,ncent,noa,nob,nvb, . mxcnfb,iconfb,dax,edb, @@ -934,7 +1018,7 @@ subroutine sfstdamat(nci,nexb,ncent,noa,nob,nvb, end subroutine sfstdamat - + subroutine sf_ptselect_uks(nexb,ncent,noa,nva,nob,nvb, . nexptb,mxcnfb,iconfb, . kconfb,dak,dax,edb,edptb, @@ -958,7 +1042,7 @@ subroutine sf_ptselect_uks(nexb,ncent,noa,nva,nob,nvb, real*8 de,pert,amat logical, allocatable :: incl_conf(:) - allocate(ptb(nexb),pt2b(nexb), + allocate(ptb(nexb),pt2b(nexb), . qj(ncent),incl_conf(nexptb),stat=ierr) if(ierr.ne.0)stop 'allocation for PT intermediates crashed' incl_conf=.false. @@ -969,8 +1053,8 @@ subroutine sf_ptselect_uks(nexb,ncent,noa,nva,nob,nvb, pt2b = 0.0d0 !$omp parallel private(i,k,io,iv,iiv,iwrk,de,j,jo,jv,jjv,jwrk,ej, -!$omp& l,qj,ptb,amat,pert) reduction (+:pt2b) -!$omp do +!$omp& l,qj,ptb,amat,pert) reduction (+:pt2b) +!$omp do c outer loop over beta S-CSF do k=1,nexptb io=kconfb(k,1) @@ -978,7 +1062,7 @@ subroutine sf_ptselect_uks(nexb,ncent,noa,nva,nob,nvb, iiv=iv-nob iwrk=(io-1)*nvb + iiv de=edptb(k) -c loop over beta P-CSF +c loop over beta P-CSF do j=1,nexb jo=iconfb(j,1) jv=iconfb(j,2) @@ -997,7 +1081,7 @@ subroutine sf_ptselect_uks(nexb,ncent,noa,nva,nob,nvb, if(pert.gt.thrp)then incl_conf(k)=.true. else -c else accumulate in the E(PT2) contribution to the alpha and beta P-CSF +c else accumulate in the E(PT2) contribution to the alpha and beta P-CSF do l=1,nexb pt2b(l)=pt2b(l) + ptb(l) enddo diff --git a/stda-rw.f b/stda-rw.f index 107815d..78254eb 100644 --- a/stda-rw.f +++ b/stda-rw.f @@ -20,15 +20,15 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ncent: number of atoms -C nmo : number of MOs +C nmo : number of MOs C nao : number of contracted AOs C xyz : array with atomic coordinates (1:3) and nuclear charge (4) in C Bohr C c : MO coefficients (nao*nmo) -C eps : MO energies (au) -C occ : occupation numbers +C eps : MO energies (au) +C occ : occupation numbers C iaoat: index array (1:nao) indicating on which atom the AO is centered -C thr : energy threshold in eV up to which energy the excited states +C thr : energy threshold in eV up to which energy the excited states C are computed = spectral range (input in eV) C thrp : threshold for perturbation selection of CSF (input) C ax : Fock exchange mixing parameter in DF used @@ -43,7 +43,7 @@ !!! used logicals from commonlogicals !!! C triplet: logical=.true. if triplet states are to be calculated C rpachk : logical=.true. if sTD-DFT is performed -C eigvec : logical=.true. print eigenvectors, ggavec is needed if eigvec and GGAs are used together +C eigvec : logical=.true. print eigenvectors, ggavec is needed if eigvec and GGAs are used together C nvec : integer, # roots for which eigenvectors are wanted C screen : prescreen in pt selection and for CSFs with small transition strengths C dokshift : shift A(ia,ia) elements if K(ia,ia) is small @@ -53,33 +53,33 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, - . ax,alphak,betaj,fthr,nvec) + . ax,alphak,betaj,fthr,nvec) use commonlogicals use commonresp use omp_lib - IMPLICIT NONE + IMPLICIT NONE c input: integer ncent,nmo,nao integer iaoat(*) - real*8 thr,thrp,ax,othr,vthr + real*8 thr,thrp,ax,othr,vthr real*8 c(*),eps(*),occ(*),xyz(4,*) logical ggavec,ex c local varibles: -c l=dipole lengths, v=dipole velocity (d/dxyz), m=angular momentum +c l=dipole lengths, v=dipole velocity (d/dxyz), m=angular momentum real*8 ::dipole(3),dipole_mo(3) integer::iwrk,jwrk real*8, allocatable ::xl(:),yl(:),zl(:) real*8, allocatable ::xv(:),yv(:),zv(:) real*8, allocatable ::xm(:),ym(:),zm(:) real*8, allocatable ::help(:),scr(:),dum(:),x(:,:) -c MOs, orbital energies and CSF printout stuff +c MOs, orbital energies and CSF printout stuff real*8, allocatable ::ca(:),epsi(:) real*8, allocatable ::umerk(:,:),umrkx(:,:),umrky(:,:),umrkz(:,:) real*8, allocatable ::rvp(:) c stuff for diag of TDA matrix (or rpa) -c critical regarding memory +c critical regarding memory integer info,lwork,liwork,il,iu,nfound real*4, allocatable ::uci (:,:) real*4, allocatable ::eci (:) @@ -92,14 +92,14 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c Linear response real*4 :: start_time, end_time, stda_time integer :: STATUS - + ccccccccccccc real*4 vu,vl integer,allocatable ::iwork(:) integer,allocatable::isuppz(:) c Löwdin MOs, repulsion terms, charges and half-transformed stuff -c critical regarding memory +c critical regarding memory real*8, allocatable ::clow(:) real*4, allocatable ::gamj(:,:) real*4, allocatable ::gamk(:,:) @@ -109,16 +109,16 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, real*4 sdot, integral !c the maximum size of the TDA expansion space - integer maxconf + integer maxconf real*8, allocatable ::ed(:),edpt(:) integer, allocatable :: iconf(:,:),kconf(:,:) - + ! check for vector printout integer, allocatable :: vecchk(:) -c intermediates +c intermediates real*8 omax,vmin,pert,de,ek,ej,ak,xc,rabx,ef,fthr,aksqrt real*8 beta1,alpha2,pp,hilf,uu,sss,rl,rv,time,coc(3) real*8 fl,fv,ec,p23,xp,umax,xvu,yvu,zvu,xmu,ymu,zmu,xlu,ylu,zlu @@ -128,17 +128,17 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, integer nci,jj,idum1,idum2,kmem,imax,jmax,lmem,ij,nroot,lin,ierr integer no,nv,n,nexpt,jhomo,nvec integer*8 imem1,imem2,imem3 -c atomic Hubbard parameters +c atomic Hubbard parameters real*8 eta(94) -c atomic masses +c atomic masses common /amass/ ams(107) real*8 ams character*79 dummy - + call cpu_time(stda_time) - -c just a printout - call header('s T D A',0) + +c just a printout + call header('s T D A',0) thr =thr /27.211385050d0 c estimate the orbital energy window which corresponds to the desired @@ -150,14 +150,14 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, moci=0 omax=-1d+42 vmin= 1d+42 - + do i=1,nmo if(occ(i).gt.1.990d0.and.eps(i).gt.omax) omax=eps(i) if(occ(i).lt.0.010d0.and.eps(i).lt.vmin)vmin=eps(i) enddo ! if eigenvectors are wanted in TM format,check now how many occupied there are in general - if(eigvec) then + if(eigvec) then jhomo=0 do i=1,nmo if(occ(i).gt.1.990d0) jhomo=jhomo+1 @@ -196,12 +196,12 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, . ca(nao*moci),epsi(moci) . ) - write(*,*)'MOs in TDA : ', moci + write(*,*)'MOs in TDA : ', moci ! make two cases: 1st one) eigenvectors are needed, 2) eigenvectors are not needed if(eigvec.or.nto) then ! we want eigenvectors to be printed out allocate(vecchk(nmo), stat=ierr) - if(ierr.ne.0)stop 'allocation failed for vecchk' + if(ierr.ne.0)stop 'allocation failed for vecchk' vecchk=0 moci=0 do i=1,nmo @@ -216,7 +216,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo ihomo=moci do i=1,nmo - if(occ(i).lt.0.010d0.and.eps(i).lt.vthr)then + if(occ(i).lt.0.010d0.and.eps(i).lt.vthr)then moci=moci+1 do j=1,nao ca(j+(moci-1)*nao)=c(j+(i-1)*nao) @@ -238,7 +238,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo ihomo=moci do i=1,nmo - if(occ(i).lt.0.010d0.and.eps(i).lt.vthr)then + if(occ(i).lt.0.010d0.and.eps(i).lt.vthr)then moci=moci+1 do j=1,nao ca(j+(moci-1)*nao)=c(j+(i-1)*nao) @@ -250,10 +250,10 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, no=ihomo nv=moci-no - write(*,*)'oMOs in TDA: ', no - write(*,*)'vMOs in TDA: ', nv - if(no.eq.0.or.nv.eq.0) then - stop 'no CSF, increase energy threshold (-e option)' + write(*,*)'oMOs in TDA: ', no + write(*,*)'vMOs in TDA: ', nv + if(no.eq.0.or.nv.eq.0) then + stop 'no CSF, increase energy threshold (-e option)' endif maxconf=no*nv @@ -265,7 +265,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, edpt=0.0d0 iconf=0 kconf=0 -c we arrange MOS according to energy from 1:HOMO to LUMO:MOCI +c we arrange MOS according to energy from 1:HOMO to LUMO:MOCI c (in TM they come in irreps) write(*,*)'sorting MOs ...' c sort for E diag @@ -284,7 +284,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, open(unit=31,file='xlint',form='unformatted',status='old') read(31) help - + ! dipole(1)=0.0 ! do k=1,moci ! Do i=1,nao @@ -295,23 +295,23 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ! enddo ! enddo ! enddo -! write(*,*)'mu_x',dipole(1),'electronic, not +! write(*,*)'mu_x',dipole(1),'electronic, not ! .shifted to the centre of mass' - + call onetri(1,help,dum,scr,ca,nao,moci) call shrink(moci,dum,xl) - + !dipole_mo(1)=0.0 !Do i=1,moci !dipole_mo(1)=dipole_mo(1)-xl(lin(i,i)) !enddo !write(*,*)dipole_mo(1) - - + + close(31,status='delete') open(unit=32,file='ylint',form='unformatted',status='old') read(32) help - + ! dipole(2)=0.0 ! do k=1,moci ! Do i=1,nao @@ -323,20 +323,20 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ! enddo ! enddo ! write(*,*)'mu_y',dipole(2) - + call onetri(1,help,dum,scr,ca,nao,moci) call shrink(moci,dum,yl) - + !dipole_mo(2)=0.0 !Do i=1,moci !dipole_mo(2)=dipole_mo(2)-yl(lin(i,i)) !enddo !write(*,*)dipole_mo(2) - + close(32,status='delete') open(unit=33,file='zlint',form='unformatted',status='old') read(33) help - + ! dipole(3)=0.0 ! do k=1,moci ! Do i=1,nao @@ -348,16 +348,16 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ! enddo ! enddo ! write(*,*)'mu_z',dipole(3) - + call onetri(1,help,dum,scr,ca,nao,moci) call shrink(moci,dum,zl) - + !dipole_mo(3)=0.0 !Do i=1,moci !dipole_mo(3)=dipole_mo(3)-zl(lin(i,i)) !enddo !write(*,*)dipole_mo(3) - + close(33,status='delete') CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC @@ -387,7 +387,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c velocity dipole c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - + open(unit=37,file='xvint',form='unformatted',status='old') read(37) help call onetri(-1,help,dum,scr,ca,nao,moci) @@ -405,18 +405,18 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, close(39,status='delete') CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -c calc S^1/2 and q(GS) +c calc S^1/2 and q(GS) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - open(unit=40,file='sint',form='unformatted',status='old') + open(unit=40,file='sint',form='unformatted',status='old') read(40) help write(*,*) 'ints done.' close(40,status='delete') - + write(*,*) 'S^1/2 ...' call makel(nao,help,x) - + call dgemm('n','n',nao,moci,nao,1.d0,X,nao,CA,nao,0.d0,SCR,nao) write(*,*) 'S^1/2 orthogonalized MO coefficients done.' @@ -437,7 +437,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ! call cpu_time(time) -c these are standard (CIS) factors for singlet +c these are standard (CIS) factors for singlet ak=2.0d0 if(triplet) ak=0.0d0 @@ -446,10 +446,10 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c close(1) c the global parameters of the method: - beta1=0.20d0 - beta2=1.830d0 - alpha1=1.420d0 - alpha2=0.480d0 + beta1=0.20d0 + beta2=1.830d0 + alpha1=1.420d0 + alpha2=0.480d0 if(betaj.lt.-99.0d0) then ! if no beta parameter was read in betaj=beta1+beta2*ax @@ -474,12 +474,12 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, xmolw=0 do i=1,ncent ii=idint(xyz(4,i)) -c ams is the atomic mass (-> mol weight for output file) +c ams is the atomic mass (-> mol weight for output file) xmolw=xmolw+ams(ii) - do j=1,i + do j=1,i jj=idint(xyz(4,j)) - xj =0.50d0*(eta(ii)+eta(jj)) * ax - xk =0.50d0*(eta(ii)+eta(jj)) + xj =0.50d0*(eta(ii)+eta(jj)) * ax + xk =0.50d0*(eta(ii)+eta(jj)) rabx=sqrt((xyz(1,i)-xyz(1,j))**2 . +(xyz(2,i)-xyz(2,j))**2 . +(xyz(3,i)-xyz(3,j))**2) @@ -489,16 +489,16 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, gamk(i,j)=gamk(j,i) enddo enddo - + ! ! write transition charges on disc ! write(*,*)'write transition charges on disc' open(unit=70,file='qii',form='unformatted',status='replace') -! open(unit=71,file='qij',form='unformatted',status='replace') - open(unit=710,file='pij',form='unformatted',status='replace') +! open(unit=71,file='qij',form='unformatted',status='replace') + open(unit=710,file='pij',form='unformatted',status='replace') allocate(q1(ncent),q2(ncent),qij(ncent,no*(no+1)/2)) - q1=0.0 + q1=0.0 q2=0.0 Do i=1, no Do j=1, i-1 @@ -507,7 +507,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, qij(1:ncent,lin(i,j))=q1(1:ncent) enddo call lo12pop(i,i,ncent,nao,iaoat,clow,q1) - write(70)q1 + write(70)q1 qij(1:ncent,lin(i,i))=q1(1:ncent) q2(1:ncent)=q2(1:ncent)+q1(1:ncent) enddo @@ -523,19 +523,19 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, deallocate(pij) close(710) open(unit=72,file='qaa',form='unformatted',status='replace') - open(unit=73,file='qab',form='unformatted',status='replace') + open(unit=73,file='qab',form='unformatted',status='replace') Do i=no+1, moci Do j=no+1, i-1 call lo12pop(i,j,ncent,nao,iaoat,clow,q1) write(73)q1 enddo call lo12pop(i,i,ncent,nao,iaoat,clow,q1) - write(72)q1 + write(72)q1 enddo close(72) close(73) open(unit=74,file='qia',form='unformatted',status='replace') - open(unit=740,file='pia',form='unformatted',status='replace') + open(unit=740,file='pia',form='unformatted',status='replace') allocate(qia(ncent,no*nv)) Do i=1, no Do j=no+1, moci @@ -543,8 +543,8 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, write(74)q1 ij=(i-1)*nv+j-no qia(1:ncent,ij)=q1(1:ncent) - enddo - enddo + enddo + enddo close(74) allocate(pia(ncent,no*nv)) pia=0.0 @@ -555,7 +555,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo close(740) deallocate(clow) - + write(*,'(/'' SCF atom population (using active MOs):'')') write(*,'(10F7.3)')q2(1:ncent)*2.0 write(*,*) @@ -563,18 +563,18 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, write(*,*) open(unit=70,file='qii',form='unformatted',status='old') - open(unit=72,file='qaa',form='unformatted',status='old') + open(unit=72,file='qaa',form='unformatted',status='old') allocate(qij(ncent,no),qab(ncent,nv)) - + Do i=1, no - read(70)qij(1:ncent,i) + read(70)qij(1:ncent,i) enddo Do i=1, nv - read(72)qab(1:ncent,i) - enddo + read(72)qab(1:ncent,i) + enddo close(70) - close(72) - + close(72) + allocate(pij(ncent,no), stat=ierr) if(ierr.ne.0)stop 'allocation failed for (ii| intermediate' pij=0.0 @@ -587,9 +587,9 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, call sgemm('t','n',nv,no,ncent,1.0,qab,ncent,pij,ncent,0.0, . uci,nv) deallocate(pij) - + allocate(q3(1:ncent)) -c determine singles which are lower than thr +c determine singles which are lower than thr k=0 j=0 do io=1,no ! occ loop @@ -597,17 +597,17 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do iv=no+1,moci ! virt loop de=epsi(iv)-epsi(io) q2(1:ncent)=qab(1:ncent,iv-no) !qaa - ej=dble(uci(iv-no,io)) + ej=dble(uci(iv-no,io)) de=de-ej i=nv*(io-1)+iv-no q2(1:ncent)=pia(1:ncent,i) - q3(1:ncent)=qia(1:ncent,i) !qia + q3(1:ncent)=qia(1:ncent,i) !qia ek=sdot(ncent,q2,1,q3,1) de=de+ak*ek ! optional: perform K(ia,ia) dependent shift if(dokshift) then - call kshift_to_ediag(de,ek) + call kshift_to_ediag(de,ek) endif c the primary ones @@ -617,7 +617,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, iconf(k,2)=iv ed(k)=de endif -c for PT +c for PT if(de.gt.thr.and.de.lt.fthr)then ! needs to be on if fthr is specified j=j+1 kconf(j,1)=io @@ -632,7 +632,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, nci=nex nex=k nexpt=j - + write(*,*) write(*,*)nex,'CSF included by energy.' write(*,*) @@ -640,32 +640,32 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c errors and warning if(nex.lt.1) stop 'no CSF, increase energy threshold (-e option)' -! if(nex.eq.maxconf) +! if(nex.eq.maxconf) ! . stop 'primary CSF space exceeded. use -e option!' -! if(nexpt.eq.maxconf) +! if(nexpt.eq.maxconf) ! . write(*,*)'CSF PT2 space exceeded. try -p option!' - + c sort for E diag - do 141 ii = 2,nex - i = ii - 1 - k = i - pp= ed(i) - do 121 j = ii, nex - if (ed(j) .ge. pp) go to 121 - k = j + do 141 ii = 2,nex + i = ii - 1 + k = i + pp= ed(i) + do 121 j = ii, nex + if (ed(j) .ge. pp) go to 121 + k = j pp=ed(j) - 121 continue - if (k .eq. i) go to 141 - ed(k) = ed(i) - ed(i) = pp - do m=1,2 - ihilf=iconf(i,m) - iconf(i,m)=iconf(k,m) - iconf(k,m)=ihilf + 121 continue + if (k .eq. i) go to 141 + ed(k) = ed(i) + ed(i) = pp + do m=1,2 + ihilf=iconf(i,m) + iconf(i,m)=iconf(k,m) + iconf(k,m)=ihilf enddo - 141 continue + 141 continue -c just printout +c just printout write(*,*)'ordered frontier orbitals' write(*,*)' eV' j=max(1,no-10) @@ -679,26 +679,26 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo open(unit=70,file='qii',form='unformatted',status='old') - open(unit=72,file='qaa',form='unformatted',status='old') + open(unit=72,file='qaa',form='unformatted',status='old') open(unit=74,file='qia',form='unformatted',status='old') allocate(qij(ncent,no),qab(ncent,nv),qia(ncent,no*nv)) - + Do i=1, no - read(70)qij(1:ncent,i) + read(70)qij(1:ncent,i) enddo Do i=1, nv - read(72)qab(1:ncent,i) + read(72)qab(1:ncent,i) enddo Do i=1, no Do j=no+1, moci ij=(i-1)*nv+j-no read(74)qia(1:ncent,ij) - enddo - enddo + enddo + enddo close(70,status='delete') close(72) - close(74) - + close(74) + write(*,*) write(*,*)' lowest CSF states' write(*,*)' eV nm excitation i->a eV' @@ -708,7 +708,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, q1(1:ncent)=qij(1:ncent,io) jii=integral(q1,q1,gamj,ncent) k=iv-no - q2(1:ncent)=qab(1:ncent,k) + q2(1:ncent)=qab(1:ncent,k) ej=integral(q1,q2,gamj,ncent) jaa=integral(q2,q2,gamj,ncent) l=nv*(io-1)+k @@ -722,7 +722,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, . 1.d+7/(ed(i)*2.19474625d+5),iconf(i,1:2), . 27.211*(epsi(iv)-epsi(io)),27.211*ej,27.211*ek,27.211*de,loc enddo - + deallocate(q1,q2,q3,qij,qab,qia) 14 format(i5,f6.2,f8.1, 5x,i4,' ->',i4,5x,'gap,J,K:',3f8.3, . 3x,'Kshft:',f8.3,2x,'locality:',f6.3,E12.5) @@ -732,7 +732,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, write(*,*) write(*,*)'selecting CSF ...' if(pt_off)then - deallocate(edpt,kconf) + deallocate(edpt,kconf) nroot=nex nci=nex write(*,*)nci,'CSF in total.' @@ -743,14 +743,14 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c nroot at this point is the number of primary CSF. The c number of roots in the range 0-thr (as desired) is not -c known but will be determined by the diag routine. +c known but will be determined by the diag routine. nroot=nex write(*,*)new,'CSF included by PT.' nci=nex+new write(*,*)nci,'CSF in total.' endif ! call cpu_time(hilf) -! write(*,*) 'time elapsed:',hilf-time +! write(*,*) 'time elapsed:',hilf-time if(rpachk) then ************************** c Obtional RPA procedure * @@ -758,87 +758,105 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, write(*,*) 'sTD-DFT procedure...' write(*,*) 'setting up A+B and A-B matrices' -c allocate A+B and A-B in packed form +c allocate A+B and A-B in packed form allocate( apb(nci*(nci+1)/2),ambsqr(nci*(nci+1)/2), - . stat=ierr ) + . stat=ierr ) if(ierr.ne.0)stop 'allocation failed for A+B or A-B' call rrpamat_rw(nci,ncent,no,nv,maxconf,iconf,ak,ax,ed,gamj . ,gamk,apb,ambsqr,moci) - + ***************************** c Linear Response functions * -***************************** +***************************** + if(optrota)then + call cpu_time(start_time) + allocate( amb(nci*(nci+1)/2), stat=ierr ) + if(ierr.ne.0)stop 'allocation failed for A-B' + open(unit=53,file='amb',form='unformatted',status='old') + read(53) amb + close(53,status='delete') + if(velo_OR==.false.)call optrot(nci,apb,amb,iconf,maxconf, + .xl,yl,zl,moci,no,nv,xm,ym,zm,xmolw) + if(velo_OR==.true.)call optrot_velo(nci,apb,amb,iconf,maxconf, + .xv,yv,zv,moci,no,nv,xm,ym,zm,xmolw) + call cpu_time(end_time) + print '("Opt. Rot. Time = ",f12.2," minutes.")' + . ,(end_time-start_time)/60.0 + print '("sTD-DFT Time = ",f12.2," minutes.")' + . ,(end_time-stda_time)/60.0 + CALL EXIT(STATUS) + endif if(resp) then if(triplet) stop 'not available' call cpu_time(start_time) - allocate( amb(nci*(nci+1)/2), stat=ierr ) - if(ierr.ne.0)stop 'allocation failed for A-B' + allocate( amb(nci*(nci+1)/2), stat=ierr ) + if(ierr.ne.0)stop 'allocation failed for A-B' open(unit=53,file='amb',form='unformatted',status='old') read(53) amb close(53,status='delete') - + call lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, . no,nv) - + call cpu_time(end_time) print '("Lresp Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 - + ! call lresp1_noinv(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . no,nv) -! +! ! call cpu_time(end_time) ! print '("Lresp Time = ",f12.2," minutes.")' ! . ,(end_time-start_time)/60.0 print '("sTD-DFT-rw Time = ",f12.2," minutes.")' . ,(end_time-stda_time)/60.0 - CALL EXIT(STATUS) + CALL EXIT(STATUS) endif - + if(aresp) then if(triplet) stop 'not available' call cpu_time(start_time) - allocate( amb(nci*(nci+1)/2), stat=ierr ) - if(ierr.ne.0)stop 'allocation failed for A-B' + allocate( amb(nci*(nci+1)/2), stat=ierr ) + if(ierr.ne.0)stop 'allocation failed for A-B' open(unit=53,file='amb',form='unformatted',status='old') read(53) amb close(53,status='delete') - + call lresp(nci,apb,ambsqr,iconf,maxconf,xl,yl,zl,moci, . no,nv) - + call cpu_time(end_time) print '("Lresp Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 - + ! call lresp_noinv(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . no,nv) ! call lresp_noinv1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . no,nv) -! +! ! call cpu_time(end_time) ! print '("Lresp Time = ",f12.2," minutes.")' ! . ,(end_time-start_time)/60.0 print '("sTD-DFT-rw Time = ",f12.2," minutes.")' . ,(end_time-stda_time)/60.0 - CALL EXIT(STATUS) + CALL EXIT(STATUS) endif - + c big arrays not needed anymore deallocate(ed) ! call prmat4(6,ambsqr,nci,0,'A-B^0.5') -************************************************************************************ +************************************************************************************ ! if eigenvectors in TM format wanted for GGAs, this is done here (not so nice) ************************************************************************************ ggavec=.false. - if (ax.eq.0.0d0.and.eigvec) then + if (ax.eq.0.0d0.and.eigvec) then open(unit=39,file='TmPvEcInFo',status='replace') ! print to temporary file - if(triplet) then + if(triplet) then write(39,*) 1 else write(39,*) 0 @@ -868,48 +886,48 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, if(info.lt.1) stop 'internal error in diag' nroot=info info=0 - - + + ! if(resp)then -! +! ! call pol_sos(nroot,nci,eci,uci,hci,xl,yl,zl,moci, ! .maxconf,iconf,ak) -! -! +! +! ! endif - + ***************************** c Lin. Response func. 2PA * -***************************** +***************************** if(TPA)then if(triplet) stop 'not available' call cpu_time(start_time) - allocate( amb(nci*(nci+1)/2), stat=ierr ) - if(ierr.ne.0)stop 'allocation failed for A-B' + allocate( amb(nci*(nci+1)/2), stat=ierr ) + if(ierr.ne.0)stop 'allocation failed for A-B' open(unit=53,file='amb',form='unformatted',status='old') read(53) amb close(53) - + call lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, . no,nv,eci,uci,hci,nroot) - + call cpu_time(end_time) print '("Lresp Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 - + ! call lresp_2PA_noinv(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . no,nv,eci,uci,hci,nroot) -! +! ! call cpu_time(end_time) ! print '("Lresp Time = ",f12.2," minutes.")' ! . ,(end_time-start_time)/60.0 - + print '("sTD-DFT-rw Time = ",f12.2," minutes.")' . ,(end_time-stda_time)/60.0 write(*,*) - !CALL EXIT([STATUS]) + !CALL EXIT([STATUS]) endif - + deallocate(apb,ambsqr) !******************************************************************************** open(unit=53,file='amb',form='unformatted',status='old') @@ -946,27 +964,27 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, write(*,'('' estimated time (min) '',f8.2)') . float(nci)**2*float(nroot)/8.d+8/60. -c if LAPACK does not work +c if LAPACK does not work c allocate(eci(nci),uci(nci,nroot), c . stat=ierr) c call sHQRII(hci,nci,nroot,eci,uci) c faster by a factor of 2-3 - lwork =26*nci + lwork =26*nci liwork=10*nci -c we allocate uci with nci (and not with nroot) as save +c we allocate uci with nci (and not with nroot) as save c choice (other values gave segfaults) nroot=min(nci,int(1.5*nroot)) allocate(eci(nci),uci(nci,nci),work(lwork), . iwork(liwork),isuppz(nci)) if(ierr.ne.0)stop 'allocation failed for TDA matrix diag' - vl=0 + vl=0 vu=thr call ssyevr('V','V','U',nci,hci,nci,vl,vu,il,iu,1.e-6, . nfound,eci,uci,nci,isuppz, . work,lwork,iwork,liwork,info) - + nroot=nfound if(nfound.lt.1) stop 'internal error in diag' @@ -985,7 +1003,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ! enddo !!!! - endif + endif write(*,'(i5,'' roots found, lowest/highest eigenvalue : '', .2F8.3,i4)')nroot,eci(1)*27.21139,eci(nroot)*27.21139,info @@ -993,7 +1011,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, allocate(umerk(14,nroot)) c contract WF with MO integrals for intensities - p23= 2.0d0/3.0d0 + p23= 2.0d0/3.0d0 aksqrt=sqrt(ak) ! if aniso, store and print the x,y,z-resolved data @@ -1003,11 +1021,11 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, umrkx=0.0d0 umrky=0.0d0 umrkz=0.0d0 - endif + endif ! if desired, print out data for exciton coupling, i.e., transition moments if(printexciton) then - open(unit=27,file='tda.exc') + open(unit=27,file='tda.exc') write(27,*)ncent,nroot write(27,'(a)',advance='yes')'$coord' ihilf=0 @@ -1028,21 +1046,21 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do i=1,nroot de=dble(eci(i)) ef=1.0d0/(de+1.0d-8) - xlu=0.0d0 - ylu=0.0d0 - zlu=0.0d0 - xmu=0.0d0 - ymu=0.0d0 - zmu=0.0d0 - xvu=0.0d0 - yvu=0.0d0 + xlu=0.0d0 + ylu=0.0d0 + zlu=0.0d0 + xmu=0.0d0 + ymu=0.0d0 + zmu=0.0d0 + xvu=0.0d0 + yvu=0.0d0 zvu=0.0d0 xms=0.0d0 yms=0.0d0 - zms=0.0d0 + zms=0.0d0 l=0 - umax=-1 - kmem=1 + umax=-1 + kmem=1 q1=0.0 c loop over CSF do j=1,nci @@ -1073,7 +1091,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ymu=ymu+ym(ij)*pp zmu=zmu+zm(ij)*pp if(.not.velcorr.and.printexciton) then - l=(io-1)*nv+(iv-no) + l=(io-1)*nv+(iv-no) q1(1:ncent)=q1(1:ncent)+qia(1:ncent,l)*real(uu) endif enddo @@ -1089,8 +1107,8 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ymu=ymu*aksqrt zmu=zmu*aksqrt -c polarizability - alp_real(1)=alp_real(1)+xlu*xlu*ef +c polarizability + alp_real(1)=alp_real(1)+xlu*xlu*ef alp_real(2)=alp_real(2)+xlu*ylu*ef alp_real(3)=alp_real(3)+ylu*ylu*ef alp_real(4)=alp_real(4)+xlu*zlu*ef @@ -1109,20 +1127,20 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do k=1,ncent write(27,'(f12.8)') real(aksqrt)*q1(k) ! factor of 2**0.5 arises from spin integration enddo -! write(27,*) +! write(27,*) endif - xp=xlu*xlu+ylu*ylu+zlu*zlu - fl=de*p23*xp - xp=xvu*xvu+yvu*yvu+zvu*zvu - fv=ef*p23*xp + xp=xlu*xlu+ylu*ylu+zlu*zlu + fl=de*p23*xp + xp=xvu*xvu+yvu*yvu+zvu*zvu + fv=ef*p23*xp xp=xlu*xmu+ylu*ymu+zlu*zmu - rl=-235.7220d0*xp -! multiplying with ak (2=singlet excitation -> restricted case) is already included in moments: see book by Harada & Nakanishi + rl=-235.7220d0*xp +! multiplying with ak (2=singlet excitation -> restricted case) is already included in moments: see book by Harada & Nakanishi ! recalculated value for Rot with recent CODATA values and changed it (used to be 235.730) ! the factor is the conversion from e*a0*mu_b in atomic units to 10^40 cgs (m_b is the Bohr magneton and a0 the Bohr radius) - xp=xvu*xmu+yvu*ymu+zvu*zmu - rv=-235.7220d0*xp*ef + xp=xvu*xmu+yvu*ymu+zvu*zmu + rv=-235.7220d0*xp*ef rvp(i)=rv c sum rule sumf=sumf+fl @@ -1198,39 +1216,39 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, endif enddo - + c sort (output) c do 140 ii = 2,nroot -c i = ii - 1 -c k = i -c pp= umerk(1,i) -c do 120 j = ii, nroot -c if (umerk(1,j) .ge. pp) go to 120 -c k = j -c pp=umerk(1,j) -c 120 continue -c if (k .eq. i) go to 140 -c umerk(1,k) = umerk(1,i) -c umerk(1,i) = pp +c i = ii - 1 +c k = i +c pp= umerk(1,i) +c do 120 j = ii, nroot +c if (umerk(1,j) .ge. pp) go to 120 +c k = j +c pp=umerk(1,j) +c 120 continue +c if (k .eq. i) go to 140 +c umerk(1,k) = umerk(1,i) +c umerk(1,i) = pp c do m=2,11 -c hilf=umerk(m,i) -c umerk(m,i)=umerk(m,k) -c umerk(m,k)=hilf +c hilf=umerk(m,i) +c umerk(m,i)=umerk(m,k) +c umerk(m,k)=hilf c enddo -c 140 continue +c 140 continue 11 format(i5,f9.3,f8.1,2f11.4,3x,3(F6.2,'(',i4,'->',i4,')')) 12 format(i5,f6.2,f8.1, 5x,i4,' ->',i4,5x,'gap,J,K:',3f8.3, . 3x,'Kshft:',f8.3) 13 format( 41x,f8.3,13x,f8.3) - + if(.not.rpachk) deallocate(hci) if(.not.velcorr) then write(*,*)'writing spectral data to tda.dat ...' - call print_tdadat(nroot,xmolw,eci,umerk(2:5,:),thr,'tda.dat') + call print_tdadat(nroot,xmolw,eci,umerk(2:5,:),thr,'tda.dat') if(aniso) then write(*,*)'writing anisotropic data ...' @@ -1242,40 +1260,40 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, call print_tdadat(nroot,xmolw,eci,umrkz(1:4,:),thr,'tdaz.dat') deallocate(umrkx,umrky,umrkz) write(*,*)'data given such that f = (f_x + f_y + f_z)/3' - write(*,*)' and R = R_x + R_y + R_z' - endif + write(*,*)' and R = R_x + R_y + R_z' + endif else - ! if velocity correction is on, get improved rotatory strengths - if (printexciton) then + ! if velocity correction is on, get improved rotatory strengths + if (printexciton) then call apbtrafoexc(nci,nroot,uci,eci,xl,yl,zl,xv,yv,zv,xm,ym,zm . ,xmolw,no,nv,coc,ncent,qia,maxconf,iconf,ak,rvp) else call apbtrafo(nci,nroot,uci,eci,xl,yl,zl,xv,yv,zv,xm,ym,zm - . ,xmolw,maxconf,iconf,ak,rvp) + . ,xmolw,maxconf,iconf,ak,rvp) endif do i=1,nroot umerk(5,i)=rvp(i) - enddo + enddo endif - if(printexciton) then + if(printexciton) then deallocate(qia) close(27) endif - -c output + +c output write(*,*) write(*,*) .'excitation energies, transition moments and TDA amplitudes' - if(velcorr) then + if(velcorr) then write(*,*) 'state eV nm fL Rv(corr)' else write(*,*) 'state eV nm fL Rv' endif do i=1,nroot ec=umerk(1,i) - write(*,11) + write(*,11) . i,ec*27.21139, . 1.d+7/(ec*2.19474625d+5),umerk(2,i),umerk(5,i) . ,umerk(6,i),idint(umerk(7,i)),idint(umerk(8,i)) @@ -1291,10 +1309,10 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, call sosor(nroot,xmolw,eci,rvp) deallocate(rvp,q1) - + ***************************** c Lin. Response func. ESA * -***************************** +***************************** if(ESA)then if(triplet) stop 'not available' ! Unrelaxed state-to-state transition dipole moments @@ -1317,9 +1335,9 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, . ,(end_time-stda_time)/60.0 endif write(*,*) - !CALL EXIT([STATUS]) + !CALL EXIT([STATUS]) endif - + ! call hyperpol_sos(nroot,nci,eci,uci,hci,xl,yl,zl,moci, ! .maxconf,iconf,no,nv) !to test the moments ! call tpa_sos(nroot,nci,eci,uci,hci,xl,yl,zl,moci, @@ -1351,13 +1369,13 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, call printvectda(rpachk,nci,nroot,uci,eci) ! TDA vectors endif endif - + if(nto)then if(rpachk)then call print_nto_rpa(uci,hci,ca,moci,nci,nroot,nao,iconf, . maxconf,no,nv) else - + call print_nto(uci,ca,moci,nci,nroot,nao,iconf,maxconf,no,nv) endif endif @@ -1366,44 +1384,44 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, if(rpachk) then write(*,*) 'sTD-DFT-rw done.' - else + else write(*,*) 'sTDA-rw done.' endif !!! old fitting part -!c used for fitting +!c used for fitting inquire(file='.REF',exist=ex) if(.not.ex) return open(unit=27,file='.REF') open(unit=28,file='.OUT') read(27,*) i,k - + if(k.eq.1)then -c take only bright ones +c take only bright ones read(27,*) hilf do j=1,nroot if(umerk(2,j).gt.0.03)then - write(28,'(4f12.6)') + write(28,'(4f12.6)') . (umerk(1,j)*27.211-hilf),hilf,umerk(1,j)*27.211,umerk(2,j) exit endif enddo elseif(k.eq.2)then -c take only very bright ones +c take only very bright ones read(27,*) hilf do j=1,nroot if(umerk(2,j).gt.0.2)then - write(28,'(4f12.6)') + write(28,'(4f12.6)') . (umerk(1,j)*27.211-hilf),hilf,umerk(1,j)*27.211,umerk(2,j) exit endif enddo elseif(k.eq.3)then -c take only very bright CD ones +c take only very bright CD ones read(27,*) hilf do j=1,nroot if(abs(umerk(5,j)).gt.100)then - write(28,'(4f12.6)') + write(28,'(4f12.6)') . (umerk(1,j)*27.211-hilf),hilf,umerk(1,j)*27.211,umerk(5,j) exit endif @@ -1412,7 +1430,7 @@ SUBROUTINE stda_rw(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do j=1,i read(27,*) hilf if(hilf.gt.0.01) - . write(28,'(4f12.6)') + . write(28,'(4f12.6)') . (umerk(1,j)*27.211-hilf),hilf,umerk(1,j)*27.211,umerk(2,j) enddo endif @@ -1439,7 +1457,7 @@ subroutine ptselect_rw(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, integer, intent(in) :: kconf(mxcnf,2),no,nv integer, intent (inout) :: iconf(mxcnf,2) integer, intent (out) :: new - real*8, intent(in) :: gamj(1:ncent,1:ncent),gamk(1:ncent,1:ncent) + real*8, intent(in) :: gamj(1:ncent,1:ncent),gamk(1:ncent,1:ncent) real*8, intent(in) :: dak,dax,thrp,edpt(mxcnf) real*8, intent(inout) :: ed(mxcnf) real*4, allocatable :: qia(:,:),pia(:,:) @@ -1452,27 +1470,27 @@ subroutine ptselect_rw(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, logical, allocatable :: incl_conf(:) ak=real(dak) ax=real(dax) - allocate(q1(ncent),q2(ncent),q3(ncent),pt(nex),pt2(nex), + allocate(q1(ncent),q2(ncent),q3(ncent),pt(nex),pt2(nex), . incl_conf(nexpt),stat=ierr) if(ierr.ne.0)stop 'allocation for PT intermediates crashed' incl_conf=.false. - - + + allocate(qia(ncent,mxcnf),qab(ncent,nv*(nv+1)/2), . pij(ncent,no*(no+1)/2),pia(ncent,mxcnf)) - + open(unit=710,file='pij',form='unformatted',status='old') - open(unit=72,file='qaa',form='unformatted',status='old') - open(unit=73,file='qab',form='unformatted',status='old') - open(unit=74,file='qia',form='unformatted',status='old') + open(unit=72,file='qaa',form='unformatted',status='old') + open(unit=73,file='qab',form='unformatted',status='old') + open(unit=74,file='qia',form='unformatted',status='old') open(unit=740,file='pia',form='unformatted',status='old') Do i=1, no Do j=1, i ij=lin(i,j) read(710)pij(1:ncent,ij) - enddo enddo - close(710) + enddo + close(710) Do i=no+1, moci k=i-no Do j=no+1, i-1 @@ -1481,7 +1499,7 @@ subroutine ptselect_rw(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, read(73)qab(1:ncent,ij) enddo ij=lin(k,k) - read(72)qab(1:ncent,ij) + read(72)qab(1:ncent,ij) enddo close(72) close(73) @@ -1490,11 +1508,11 @@ subroutine ptselect_rw(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, ij=(i-1)*nv+j-no read(74)qia(1:ncent,ij) read(740)pia(1:ncent,ij) - enddo - enddo - close(74) + enddo + enddo + close(74) close(740) - + new=0 pt2=0.0d0 @@ -1502,7 +1520,7 @@ subroutine ptselect_rw(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, ! loop over secondary/neglected CSF, this is done omp-parallel !$omp parallel private(i,k,io,iv,iiv,iwrk,q1,de,j,jo,jv,jjv,jwrk,ek,ej, !$omp& q2,q3,pt,amat,pert) reduction (+:pt2) -!$omp do +!$omp do do k=1,nexpt io=kconf(k,1) iv=kconf(k,2) @@ -1510,7 +1528,7 @@ subroutine ptselect_rw(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, iwrk=(io-1)*nv + iiv q1(1:ncent)=pia(1:ncent,iwrk) de=edpt(k) -c loop over primary CSF +c loop over primary CSF do j=1,nex jo=iconf(j,1) jv=iconf(j,2) @@ -1523,16 +1541,16 @@ subroutine ptselect_rw(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, q3(1:ncent)=qab(1:ncent,lin(iiv,jjv)) ej=sdot(ncent,q2,1,q3,1) pt(j)=0.0d0 - amat=ak*ek-ej + amat=ak*ek-ej pt(j)=amat**2/(de-ed(j)+1.d-10) enddo pert=sum(pt(1:nex)) -c if sum > threshold include it +c if sum > threshold include it if(pert.gt.thrp)then incl_conf(k)=.true. else do i=1,nex - pt2(i)=pt2(i)+pt(i) + pt2(i)=pt2(i)+pt(i) enddo endif enddo @@ -1542,19 +1560,19 @@ subroutine ptselect_rw(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, do i=1,nexpt if(.not.incl_conf(i)) cycle io=kconf(i,1) - iv=kconf(i,2) - de=edpt(i) + iv=kconf(i,2) + de=edpt(i) new=new+1 - iconf(nex+new,1)=io - iconf(nex+new,2)=iv + iconf(nex+new,1)=io + iconf(nex+new,2)=iv ed (nex+new )=de enddo pert=0.0d0 do i=1,nex - pert=pert+pt2(i) - ed(i)=ed(i)-pt2(i) + pert=pert+pt2(i) + ed(i)=ed(i)-pt2(i) enddo - + write(*,'('' average/max PT2 energy lowering (eV):'',2F10.3)') . 27.21139*pert/float(nex),maxval(pt2)*27.21139 @@ -1575,7 +1593,7 @@ subroutine rrpamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, use omp_lib implicit none integer, intent(in) :: nci,ncent,no,nv,mxcnf,iconf(mxcnf,2) - real*8, intent(in) :: gamj(1:ncent,1:ncent),gamk(1:ncent,1:ncent) + real*8, intent(in) :: gamj(1:ncent,1:ncent),gamk(1:ncent,1:ncent) real*4, allocatable :: qia(:,:),pia(:,:) real*4, allocatable :: qab(:,:),qij(:,:),pij(:,:) real*8, intent(in) :: dak,dax,ed(mxcnf) @@ -1592,19 +1610,19 @@ subroutine rrpamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, ambsqr=0.0e0 allocate(qab(ncent,nv*(nv+1)/2), . pij(ncent,no*(no+1)/2)) - + open(unit=710,file='pij',form='unformatted',status='old') - open(unit=72,file='qaa',form='unformatted',status='old') - open(unit=73,file='qab',form='unformatted',status='old') - open(unit=74,file='qia',form='unformatted',status='old') + open(unit=72,file='qaa',form='unformatted',status='old') + open(unit=73,file='qab',form='unformatted',status='old') + open(unit=74,file='qia',form='unformatted',status='old') open(unit=740,file='pia',form='unformatted',status='old') Do i=1, no Do j=1, i ij=lin(i,j) read(710)pij(1:ncent,ij) - enddo enddo - close(710,status='delete') + enddo + close(710,status='delete') Do i=no+1, moci k=i-no Do j=no+1, i-1 @@ -1613,7 +1631,7 @@ subroutine rrpamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, read(73)qab(1:ncent,ij) enddo ij=lin(k,k) - read(72)qab(1:ncent,ij) + read(72)qab(1:ncent,ij) enddo close(72,status='delete') close(73,status='delete') @@ -1638,22 +1656,22 @@ subroutine rrpamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, enddo ! j enddo ! i !$omp end do -!$omp end parallel +!$omp end parallel endif deallocate(pij,qab) allocate(qia(ncent,mxcnf), - . pia(ncent,mxcnf)) + . pia(ncent,mxcnf)) Do i=1, no Do j=no+1, moci ij=(i-1)*nv+j-no read(74)qia(1:ncent,ij) read(740)pia(1:ncent,ij) - enddo - enddo - close(74,status='delete') + enddo + enddo + close(74,status='delete') close(740,status='delete') - - + + ! if ax=0, A-B is diagonal and its off-diagonal elements do not need to be calculated if(abs(dax).lt.1.0d-6) then ij=0 @@ -1676,7 +1694,7 @@ subroutine rrpamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, ek=sdot(ncent,q1,1,q2,1)! ek = (ia|jb) apb(ij)=2.0*ak*ek ambsqr(ij)=0.0 - enddo ! j + enddo ! j de=real(ed(i)) ij=lin(i,i) q2(1:ncent)=qia(1:ncent,iwrk) @@ -1685,8 +1703,8 @@ subroutine rrpamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, apb(ij)=de+ak*ek ! diagonal element of A+B enddo ! i !$omp end do -!$omp end parallel - else +!$omp end parallel + else ij=0 ! for now ambsqr=A+B and apb=A-B, since we need to take the sqrt of A-B (but want to save memory) !$omp parallel private(ij,i,j,io,iv,jo,jv,iiv,iwrk,jjv, @@ -1729,28 +1747,28 @@ subroutine rrpamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, ! call prmat4(6,apb,nci,0,'A-B') ! call prmat4(6,ambsqr,nci,0,'A+B') - if(aresp.or.resp) then + if(aresp.or.resp.or.optrota) then write(*,*) ' calculating (A-B)^0.5 not necessary...' - open(unit=53,file='amb',form='unformatted',status='replace') + open(unit=53,file='amb',form='unformatted',status='replace') write(53) apb close(53) apb=ambsqr else - open(unit=52,file='apbmat',form='unformatted',status='replace') + open(unit=52,file='apbmat',form='unformatted',status='replace') write(52) ambsqr - open(unit=53,file='amb',form='unformatted',status='replace') + open(unit=53,file='amb',form='unformatted',status='replace') write(53) apb write(*,*) ' calculating (A-B)^0.5 ...' write(*,'('' estimated time (min) '',f8.2)') . float(nci)**2*float(nci)/4.d+8/60. - call smatpow(nci,apb) ! calculate sqrt(a-b), termed ambsqr + call smatpow(nci,apb) ! calculate sqrt(a-b), termed ambsqr ambsqr=apb rewind(52) read(52) apb close(52,status='delete') - close(53) + close(53) endif @@ -1768,7 +1786,7 @@ subroutine rtdamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax use omp_lib implicit none integer, intent(in) :: nci,ncent,no,nv,mxcnf,iconf(mxcnf,2) - real*8, intent(in) :: gamj(1:ncent,1:ncent),gamk(1:ncent,1:ncent) + real*8, intent(in) :: gamj(1:ncent,1:ncent),gamk(1:ncent,1:ncent) real*4, allocatable :: qia(:,:),pia(:,:) real*4, allocatable :: qab(:,:),qij(:,:),pij(:,:) real*8, intent(in) :: dak,dax,ed(mxcnf) @@ -1780,22 +1798,22 @@ subroutine rtdamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax if(ierr.ne.0)stop 'allocation for qkj crashed' ak=real(dak) ax=real(dax) - + allocate(qia(ncent,mxcnf),qab(ncent,nv*(nv+1)/2), . pij(ncent,no*(no+1)/2),pia(ncent,mxcnf)) - + open(unit=710,file='pij',form='unformatted',status='old') - open(unit=72,file='qaa',form='unformatted',status='old') - open(unit=73,file='qab',form='unformatted',status='old') - open(unit=74,file='qia',form='unformatted',status='old') + open(unit=72,file='qaa',form='unformatted',status='old') + open(unit=73,file='qab',form='unformatted',status='old') + open(unit=74,file='qia',form='unformatted',status='old') open(unit=740,file='pia',form='unformatted',status='old') Do i=1, no Do j=1, i ij=lin(i,j) read(710)pij(1:ncent,ij) - enddo enddo - close(710,status='delete') + enddo + close(710,status='delete') Do i=no+1, moci k=i-no Do j=no+1, i-1 @@ -1804,7 +1822,7 @@ subroutine rtdamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax read(73)qab(1:ncent,ij) enddo ij=lin(k,k) - read(72)qab(1:ncent,ij) + read(72)qab(1:ncent,ij) enddo close(72,status='delete') close(73,status='delete') @@ -1813,11 +1831,11 @@ subroutine rtdamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax ij=(i-1)*nv+j-no read(74)qia(1:ncent,ij) read(740)pia(1:ncent,ij) - enddo - enddo - close(74,status='delete') + enddo + enddo + close(74,status='delete') close(740,status='delete') - + ! calculate CIS matrix A hci=0.0e0 !$omp parallel private(i,j,io,iv,jo,jv,iiv,iwrk,jjv,jwrk,q1, @@ -1842,16 +1860,16 @@ subroutine rtdamat_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax hci(j,i)=ak*ek-ej hci(i,j)=hci(j,i) enddo - hci(i,i)=real(ed(i)) + hci(i,i)=real(ed(i)) enddo !$omp end do !$omp end parallel - deallocate(qia,qab,pij,pia) + deallocate(qia,qab,pij,pia) deallocate(q1,q2,q3) return end subroutine rtdamat_rw -*********************************************************************** - +*********************************************************************** + real*4 function integral(q_ij,q_ab,gamma,ncent) use omp_lib implicit none @@ -1864,7 +1882,7 @@ real*4 function integral(q_ij,q_ab,gamma,ncent) integral = sdot(ncent,q_ij,1,intermediate,1) return end - + *********************************************************************** * set up 0.5*B (packed form) in RKS case *********************************************************************** @@ -1873,7 +1891,7 @@ subroutine rtdacorr_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax use omp_lib implicit none integer, intent(in) :: nci,ncent,no,nv,mxcnf,iconf(mxcnf,2) - real*8, intent(in) :: gamj(1:ncent,1:ncent),gamk(1:ncent,1:ncent) + real*8, intent(in) :: gamj(1:ncent,1:ncent),gamk(1:ncent,1:ncent) real*4, allocatable :: qia(:,:),pia(:,:) real*4, allocatable :: qab(:,:),qij(:,:),pij(:,:) real*8, intent(in) :: dak,dax,ed(mxcnf) @@ -1887,22 +1905,22 @@ subroutine rtdacorr_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax if(ierr.ne.0)stop 'allocation for qkj/bmat crashed' ak=real(dak) ax=real(dax) - + allocate(qia(ncent,mxcnf),qab(ncent,nv*(nv+1)/2), . pij(ncent,no*(no+1)/2),pia(ncent,mxcnf)) - + open(unit=710,file='pij',form='unformatted',status='old') - open(unit=72,file='qaa',form='unformatted',status='old') - open(unit=73,file='qab',form='unformatted',status='old') - open(unit=74,file='qia',form='unformatted',status='old') + open(unit=72,file='qaa',form='unformatted',status='old') + open(unit=73,file='qab',form='unformatted',status='old') + open(unit=74,file='qia',form='unformatted',status='old') open(unit=740,file='pia',form='unformatted',status='old') Do i=1, no Do j=1, i ij=lin8(i,j) read(710)pij(1:ncent,ij) - enddo enddo - close(710) + enddo + close(710) Do i=no+1, moci k=i-no Do j=no+1, i-1 @@ -1911,7 +1929,7 @@ subroutine rtdacorr_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax read(73)qab(1:ncent,ij) enddo ij=lin8(k,k) - read(72)qab(1:ncent,ij) + read(72)qab(1:ncent,ij) enddo close(72) close(73) @@ -1920,15 +1938,15 @@ subroutine rtdacorr_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax ij=(i-1)*nv+j-no read(74)qia(1:ncent,ij) read(740)pia(1:ncent,ij) - enddo - enddo - close(74) + enddo + enddo + close(74) close(740) - - -! calculate 0.5*B + + +! calculate 0.5*B bmat=0.0e0 - fact=0.50d0 ! this is the scaling of the B-contribution + fact=0.50d0 ! this is the scaling of the B-contribution open(unit=52,file='bmat',form='unformatted',status='replace') ij=0 !$omp parallel private(ij,i,j,io,iv,jo,jv,iiv,iwrk,jjv, @@ -1962,10 +1980,10 @@ subroutine rtdacorr_rw(nci,ncent,no,nv,mxcnf,iconf,dak,dax bmat(ij)=fact*(ak*ek-ax*ek) ! diagonal element of 0.5*B enddo !$omp end do -!$omp end parallel +!$omp end parallel write(52)bmat close(52) - deallocate(qia,qab,pij,pia) + deallocate(qia,qab,pij,pia) deallocate(bmat,q1,q2,q3) return diff --git a/stda.f b/stda.f index b02268d..429e342 100644 --- a/stda.f +++ b/stda.f @@ -17,15 +17,15 @@ ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ncent: number of atoms -C nmo : number of MOs +C nmo : number of MOs C nao : number of contracted AOs C xyz : array with atomic coordinates (1:3) and nuclear charge (4) in C Bohr C c : MO coefficients (nao*nmo) -C eps : MO energies (au) -C occ : occupation numbers +C eps : MO energies (au) +C occ : occupation numbers C iaoat: index array (1:nao) indicating on which atom the AO is centered -C thr : energy threshold in eV up to which energy the excited states +C thr : energy threshold in eV up to which energy the excited states C are computed = spectral range (input in eV) C thrp : threshold for perturbation selection of CSF (input) C ax : Fock exchange mixing parameter in DF used @@ -40,7 +40,7 @@ !!! used logicals from commonlogicals !!! C triplet: logical=.true. if triplet states are to be calculated C rpachk : logical=.true. if sTD-DFT is performed -C eigvec : logical=.true. print eigenvectors, ggavec is needed if eigvec and GGAs are used together +C eigvec : logical=.true. print eigenvectors, ggavec is needed if eigvec and GGAs are used together C nvec : integer, # roots for which eigenvectors are wanted C screen : prescreen in pt selection and for CSFs with small transition strengths C dokshift : shift A(ia,ia) elements if K(ia,ia) is small @@ -50,33 +50,33 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, - . ax,alphak,betaj,fthr,nvec) + . ax,alphak,betaj,fthr,nvec) use commonlogicals use commonresp use omp_lib - IMPLICIT NONE + IMPLICIT NONE c input: integer ncent,nmo,nao integer iaoat(*) - real*8 thr,thrp,ax,othr,vthr + real*8 thr,thrp,ax,othr,vthr real*8 c(*),eps(*),occ(*),xyz(4,*) logical ggavec,ex c local varibles: -c l=dipole lengths, v=dipole velocity (d/dxyz), m=angular momentum +c l=dipole lengths, v=dipole velocity (d/dxyz), m=angular momentum real*8 ::dipole(3),dipole_mo(3) integer::iwrk,jwrk real*8, allocatable ::xl(:),yl(:),zl(:) real*8, allocatable ::xv(:),yv(:),zv(:) real*8, allocatable ::xm(:),ym(:),zm(:) real*8, allocatable ::help(:),scr(:),dum(:),x(:,:) -c MOs, orbital energies and CSF printout stuff +c MOs, orbital energies and CSF printout stuff real*8, allocatable ::ca(:),epsi(:) real*8, allocatable ::umerk(:,:),umrkx(:,:),umrky(:,:),umrkz(:,:) real*8, allocatable ::rvp(:) c stuff for diag of TDA matrix (or rpa) -c critical regarding memory +c critical regarding memory integer info,lwork,liwork,il,iu,nfound real*4, allocatable ::uci (:,:) real*4, allocatable ::eci (:) @@ -89,14 +89,14 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c Linear response real*4 :: start_time, end_time, stda_time integer :: STATUS - + ccccccccccccc real*4 vu,vl integer,allocatable ::iwork(:) integer,allocatable::isuppz(:) c Löwdin MOs, repulsion terms, charges and half-transformed stuff -c critical regarding memory +c critical regarding memory real*8, allocatable ::clow(:) real*4, allocatable ::gamj(:,:) real*4, allocatable ::gamk(:,:) @@ -106,16 +106,16 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, real*4 sdot !c the maximum size of the TDA expansion space - integer maxconf + integer maxconf real*8, allocatable ::ed(:),edpt(:) integer, allocatable :: iconf(:,:),kconf(:,:) - + ! check for vector printout integer, allocatable :: vecchk(:) -c intermediates +c intermediates real*8 omax,vmin,pert,de,ek,ej,ak,xc,rabx,ef,fthr,aksqrt real*8 beta1,alpha2,pp,hilf,uu,sss,rl,rv,time,coc(3) real*8 fl,fv,ec,p23,xp,umax,xvu,yvu,zvu,xmu,ymu,zmu,xlu,ylu,zlu @@ -125,17 +125,17 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, integer nci,jj,idum1,idum2,kmem,imax,jmax,lmem,ij,nroot,lin,ierr integer no,nv,n,nexpt,jhomo,nvec integer*8 imem1,imem2,imem3 -c atomic Hubbard parameters +c atomic Hubbard parameters real*8 eta(94) -c atomic masses +c atomic masses common /amass/ ams(107) real*8 ams character*79 dummy - + call cpu_time(stda_time) - -c just a printout - call header('s T D A',0) + +c just a printout + call header('s T D A',0) thr =thr /27.211385050d0 c estimate the orbital energy window which corresponds to the desired @@ -147,14 +147,14 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, moci=0 omax=-1d+42 vmin= 1d+42 - + do i=1,nmo if(occ(i).gt.1.990d0.and.eps(i).gt.omax) omax=eps(i) if(occ(i).lt.0.010d0.and.eps(i).lt.vmin)vmin=eps(i) enddo ! if eigenvectors are wanted in TM format,check now how many occupied there are in general - if(eigvec) then + if(eigvec) then jhomo=0 do i=1,nmo if(occ(i).gt.1.990d0) jhomo=jhomo+1 @@ -193,12 +193,12 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, . ca(nao*moci),epsi(moci) . ) - write(*,*)'MOs in TDA : ', moci + write(*,*)'MOs in TDA : ', moci ! make two cases: 1st one) eigenvectors are needed, 2) eigenvectors are not needed if(eigvec.or.nto) then ! we want eigenvectors to be printed out allocate(vecchk(nmo), stat=ierr) - if(ierr.ne.0)stop 'allocation failed for vecchk' + if(ierr.ne.0)stop 'allocation failed for vecchk' vecchk=0 moci=0 do i=1,nmo @@ -213,7 +213,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo ihomo=moci do i=1,nmo - if(occ(i).lt.0.010d0.and.eps(i).lt.vthr)then + if(occ(i).lt.0.010d0.and.eps(i).lt.vthr)then moci=moci+1 do j=1,nao ca(j+(moci-1)*nao)=c(j+(i-1)*nao) @@ -235,7 +235,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo ihomo=moci do i=1,nmo - if(occ(i).lt.0.010d0.and.eps(i).lt.vthr)then + if(occ(i).lt.0.010d0.and.eps(i).lt.vthr)then moci=moci+1 do j=1,nao ca(j+(moci-1)*nao)=c(j+(i-1)*nao) @@ -247,10 +247,10 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, no=ihomo nv=moci-no - write(*,*)'oMOs in TDA: ', no - write(*,*)'vMOs in TDA: ', nv - if(no.eq.0.or.nv.eq.0) then - stop 'no CSF, increase energy threshold (-e option)' + write(*,*)'oMOs in TDA: ', no + write(*,*)'vMOs in TDA: ', nv + if(no.eq.0.or.nv.eq.0) then + stop 'no CSF, increase energy threshold (-e option)' endif maxconf=no*nv @@ -262,7 +262,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, edpt=0.0d0 iconf=0 kconf=0 -c we arrange MOS according to energy from 1:HOMO to LUMO:MOCI +c we arrange MOS according to energy from 1:HOMO to LUMO:MOCI c (in TM they come in irreps) write(*,*)'sorting MOs ...' c sort for E diag @@ -322,7 +322,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c velocity dipole c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - + open(unit=37,file='xvint',form='unformatted',status='old') read(37) help call onetri(-1,help,dum,scr,ca,nao,moci) @@ -340,18 +340,18 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, close(39,status='delete') CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -c calc S^1/2 and q(GS) +c calc S^1/2 and q(GS) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - open(unit=40,file='sint',form='unformatted',status='old') + open(unit=40,file='sint',form='unformatted',status='old') read(40) help write(*,*) 'ints done.' close(40,status='delete') - + write(*,*) 'S^1/2 ...' call makel(nao,help,x) - + call dgemm('n','n',nao,moci,nao,1.d0,X,nao,CA,nao,0.d0,SCR,nao) write(*,*) 'S^1/2 orthogonalized MO coefficients done.' @@ -378,7 +378,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, if(ierr.ne.0) stop 'error in diag. J charges allocation' qij=0.0 qab=0.0 - + do i=1,ihomo call lo12pop(i,i,ncent,nao,iaoat,clow,q2) q1(1:ncent)=q1(1:ncent)+q2(1:ncent) @@ -395,7 +395,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, qab(1:ncent,j)=q2(1:ncent) enddo -c these are standard (CIS) factors for singlet +c these are standard (CIS) factors for singlet ak=2.0d0 if(triplet) ak=0.0d0 @@ -404,10 +404,10 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c close(1) c the global parameters of the method: - beta1=0.20d0 - beta2=1.830d0 - alpha1=1.420d0 - alpha2=0.480d0 + beta1=0.20d0 + beta2=1.830d0 + alpha1=1.420d0 + alpha2=0.480d0 if(betaj.lt.-99.0d0) then ! if no beta parameter was read in betaj=beta1+beta2*ax @@ -432,12 +432,12 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, xmolw=0 do i=1,ncent ii=idint(xyz(4,i)) -c ams is the atomic mass (-> mol weight for output file) +c ams is the atomic mass (-> mol weight for output file) xmolw=xmolw+ams(ii) - do j=1,i + do j=1,i jj=idint(xyz(4,j)) - xj =0.50d0*(eta(ii)+eta(jj)) * ax - xk =0.50d0*(eta(ii)+eta(jj)) + xj =0.50d0*(eta(ii)+eta(jj)) * ax + xk =0.50d0*(eta(ii)+eta(jj)) rabx=sqrt((xyz(1,i)-xyz(1,j))**2 . +(xyz(2,i)-xyz(2,j))**2 . +(xyz(3,i)-xyz(3,j))**2) @@ -466,8 +466,8 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, -! compute intermediates q's refer to charges, p's to gam*q (i.e., contracted) - +! compute intermediates q's refer to charges, p's to gam*q (i.e., contracted) + ! Coulomb type terms: calc qij*gam^J allocate(pij(ncent,no), stat=ierr) if(ierr.ne.0)stop 'allocation failed for (ii| intermediate' @@ -488,13 +488,13 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, if(ierr.ne.0)stop 'allocation failed for intermediates' qia=0.0 pia=0.0 -! K-type terms +! K-type terms k=0 !$omp parallel private(i,j,k,q2) !$omp do do i=1,no - do j=no+1,moci -! k=k+1 + do j=no+1,moci +! k=k+1 k=(i-1)*nv+j-no call lo12pop(i,j,ncent,nao,iaoat,clow,q2) qia(1:ncent,k)=q2(1:ncent) @@ -502,12 +502,12 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo !$omp end do !$omp end parallel - pia=0.0 + pia=0.0 call ssymm('l','l',ncent,nex,1.0,gamk,ncent,qia,ncent,0.0 . ,pia,ncent) deallocate(gamk) -c determine singles which are lower than thr +c determine singles which are lower than thr k=0 j=0 l=0 @@ -521,13 +521,13 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, de=de-ej ek=0.0d0 i=nv*(io-1)+l - q1(1:ncent)=pia(1:ncent,i) + q1(1:ncent)=pia(1:ncent,i) ek=sdot(ncent,q1,1,qia(1,i),1) de=de+ak*ek ! optional: perform K(ia,ia) dependent shift if(dokshift) then - call kshift_to_ediag(de,ek) + call kshift_to_ediag(de,ek) endif c the primary ones @@ -537,7 +537,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, iconf(k,2)=iv ed(k)=de endif -c for PT +c for PT if(de.gt.thr.and.de.lt.fthr)then ! needs to be on if fthr is specified j=j+1 kconf(j,1)=io @@ -551,7 +551,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, nci=nex nex=k nexpt=j - + write(*,*) write(*,*)nex,'CSF included by energy.' write(*,*) @@ -559,32 +559,32 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c errors and warning if(nex.lt.1) stop 'no CSF, increase energy threshold (-e option)' -! if(nex.eq.maxconf) +! if(nex.eq.maxconf) ! . stop 'primary CSF space exceeded. use -e option!' -! if(nexpt.eq.maxconf) +! if(nexpt.eq.maxconf) ! . write(*,*)'CSF PT2 space exceeded. try -p option!' - + c sort for E diag - do 141 ii = 2,nex - i = ii - 1 - k = i - pp= ed(i) - do 121 j = ii, nex - if (ed(j) .ge. pp) go to 121 - k = j + do 141 ii = 2,nex + i = ii - 1 + k = i + pp= ed(i) + do 121 j = ii, nex + if (ed(j) .ge. pp) go to 121 + k = j pp=ed(j) - 121 continue - if (k .eq. i) go to 141 - ed(k) = ed(i) - ed(i) = pp - do m=1,2 - ihilf=iconf(i,m) - iconf(i,m)=iconf(k,m) - iconf(k,m)=ihilf + 121 continue + if (k .eq. i) go to 141 + ed(k) = ed(i) + ed(i) = pp + do m=1,2 + ihilf=iconf(i,m) + iconf(i,m)=iconf(k,m) + iconf(k,m)=ihilf enddo - 141 continue + 141 continue -c just printout +c just printout write(*,*)'ordered frontier orbitals' write(*,*)' eV # centers' j=max(1,no-10) @@ -593,7 +593,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do k=1,ncent xc=xc+qij(k,i)**2 enddo - write(*,'(i4,F10.3,F8.1)') i,epsi(i)*27.21139,1./(xc+1.d-8) + write(*,'(i4,F10.3,F8.1)') i,epsi(i)*27.21139,1./(xc+1.d-8) enddo write(*,*) j=min(moci,no+11) @@ -602,7 +602,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do k=1,ncent xc=xc+qab(k,i-no)**2 enddo - write(*,'(i4,F10.3,F8.1)') i,epsi(i)*27.21139,1./(xc+1.d-8) + write(*,'(i4,F10.3,F8.1)') i,epsi(i)*27.21139,1./(xc+1.d-8) enddo @@ -616,13 +616,13 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, call ssymv('l',ncent,1.0e0,gamj,ncent,q1,1,0.0,q2,1) jii=sdot(ncent,q1,1,q2,1) k=iv-no - q1(1:ncent)=qab(1:ncent,k) + q1(1:ncent)=qab(1:ncent,k) ej=sdot(ncent,q1,1,q2,1) call ssymv('l',ncent,1.0e0,gamj,ncent,q1,1,0.0,q2,1) jaa=sdot(ncent,q1,1,q2,1) l=nv*(io-1)+k q1(1:ncent)=pia(1:ncent,l) - q2(1:ncent)=qia(1:ncent,l) + q2(1:ncent)=qia(1:ncent,l) ek=sdot(ncent,q1,1,q2,1) ! de is now the Kia shift de=0 @@ -652,14 +652,14 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo enddo !$omp end do -!$omp end parallel +!$omp end parallel allocate(pij(ncent,ihilf), stat=ierr) if(ierr.ne.0)stop 'allocation failed for (ij| intermediate' pij=0.0 call ssymm('l','l',ncent,ihilf,1.0,gamj,ncent,qij,ncent,0.0 . ,pij,ncent) - + q2=0.0 ihilf=nv*(nv+1)/2 allocate(qab(ncent,ihilf),stat=ierr) @@ -677,8 +677,8 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, enddo enddo !$omp end do -!$omp end parallel - deallocate(clow) ! orthogonalized MO coefficients not needed any more +!$omp end parallel + deallocate(clow) ! orthogonalized MO coefficients not needed any more ! call cpu_time(hilf) ! write(*,*) 'time elapsed:',hilf-time @@ -692,14 +692,14 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c nroot at this point is the number of primary CSF. The c number of roots in the range 0-thr (as desired) is not -c known but will be determined by the diag routine. +c known but will be determined by the diag routine. nroot=nex write(*,*)new,'CSF included by PT.' nci=nex+new write(*,*)nci,'CSF in total.' ! call cpu_time(hilf) -! write(*,*) 'time elapsed:',hilf-time +! write(*,*) 'time elapsed:',hilf-time if(rpachk) then ************************** c Obtional RPA procedure * @@ -707,88 +707,106 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, write(*,*) 'sTD-DFT procedure...' write(*,*) 'setting up A+B and A-B matrices' -c allocate A+B and A-B in packed form +c allocate A+B and A-B in packed form allocate( apb(nci*(nci+1)/2),ambsqr(nci*(nci+1)/2), - . stat=ierr ) + . stat=ierr ) if(ierr.ne.0)stop 'allocation failed for A+B or A-B' call rrpamat(nci,ncent,no,nv,maxconf,iconf,ak,ax,ed,pia . ,qia,pij,qab,apb,ambsqr) - + ***************************** c Linear Response functions * -***************************** +***************************** + if(optrota)then + call cpu_time(start_time) + allocate( amb(nci*(nci+1)/2), stat=ierr ) + if(ierr.ne.0)stop 'allocation failed for A-B' + open(unit=53,file='amb',form='unformatted',status='old') + read(53) amb + close(53,status='delete') + if(velo_OR==.false.)call optrot(nci,apb,amb,iconf,maxconf, + .xl,yl,zl,moci,no,nv,xm,ym,zm,xmolw) + if(velo_OR==.true.)call optrot_velo(nci,apb,amb,iconf,maxconf, + .xv,yv,zv,moci,no,nv,xm,ym,zm,xmolw) + call cpu_time(end_time) + print '("Opt. Rot. Time = ",f12.2," minutes.")' + . ,(end_time-start_time)/60.0 + print '("sTD-DFT Time = ",f12.2," minutes.")' + . ,(end_time-stda_time)/60.0 + CALL EXIT(STATUS) + endif if(resp) then if(triplet) stop 'not available' call cpu_time(start_time) - allocate( amb(nci*(nci+1)/2), stat=ierr ) - if(ierr.ne.0)stop 'allocation failed for A-B' + allocate( amb(nci*(nci+1)/2), stat=ierr ) + if(ierr.ne.0)stop 'allocation failed for A-B' open(unit=53,file='amb',form='unformatted',status='old') read(53) amb close(53,status='delete') - + call lresp1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, . no,nv) - + call cpu_time(end_time) print '("Lresp Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 - + ! call lresp1_noinv(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . no,nv) -! +! ! call cpu_time(end_time) ! print '("Lresp Time = ",f12.2," minutes.")' ! . ,(end_time-start_time)/60.0 print '("sTD-DFT Time = ",f12.2," minutes.")' . ,(end_time-stda_time)/60.0 - CALL EXIT(STATUS) + CALL EXIT(STATUS) endif - + if(aresp) then if(triplet) stop 'not available' call cpu_time(start_time) - allocate( amb(nci*(nci+1)/2), stat=ierr ) - if(ierr.ne.0)stop 'allocation failed for A-B' + allocate( amb(nci*(nci+1)/2), stat=ierr ) + if(ierr.ne.0)stop 'allocation failed for A-B' open(unit=53,file='amb',form='unformatted',status='old') read(53) amb close(53,status='delete') - + call lresp(nci,apb,ambsqr,iconf,maxconf,xl,yl,zl,moci, . no,nv) - + call cpu_time(end_time) print '("Lresp Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 - + ! call lresp_noinv(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . no,nv) ! call lresp_noinv1(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . no,nv) -! +! ! call cpu_time(end_time) ! print '("Lresp Time = ",f12.2," minutes.")' ! . ,(end_time-start_time)/60.0 print '("sTD-DFT Time = ",f12.2," minutes.")' . ,(end_time-stda_time)/60.0 - CALL EXIT(STATUS) + CALL EXIT(STATUS) endif - + c big arrays not needed anymore deallocate(pij,qab,pia,ed) if(.not.printexciton)deallocate(qia) ! call prmat4(6,ambsqr,nci,0,'A-B^0.5') -************************************************************************************ +************************************************************************************ ! if eigenvectors in TM format wanted for GGAs, this is done here (not so nice) ************************************************************************************ ggavec=.false. - if (ax.eq.0.0d0.and.eigvec) then + if (ax.eq.0.0d0.and.eigvec) then open(unit=39,file='TmPvEcInFo',status='replace') ! print to temporary file - if(triplet) then + if(triplet) then write(39,*) 1 else write(39,*) 0 @@ -818,48 +836,48 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, if(info.lt.1) stop 'internal error in diag' nroot=info info=0 - - + + ! if(resp)then -! +! ! call pol_sos(nroot,nci,eci,uci,hci,xl,yl,zl,moci, ! .maxconf,iconf,ak) -! -! +! +! ! endif - + ***************************** c Lin. Response func. 2PA * -***************************** +***************************** if(TPA)then if(triplet) stop 'not available' call cpu_time(start_time) - allocate( amb(nci*(nci+1)/2), stat=ierr ) - if(ierr.ne.0)stop 'allocation failed for A-B' + allocate( amb(nci*(nci+1)/2), stat=ierr ) + if(ierr.ne.0)stop 'allocation failed for A-B' open(unit=53,file='amb',form='unformatted',status='old') read(53) amb close(53) - + call lresp_2PA(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, . no,nv,eci,uci,hci,nroot) - + call cpu_time(end_time) print '("Lresp Time = ",f12.2," minutes.")' . ,(end_time-start_time)/60.0 - + ! call lresp_2PA_noinv(nci,apb,amb,iconf,maxconf,xl,yl,zl,moci, ! . no,nv,eci,uci,hci,nroot) -! +! ! call cpu_time(end_time) ! print '("Lresp Time = ",f12.2," minutes.")' ! . ,(end_time-start_time)/60.0 - + print '("sTD-DFT Time = ",f12.2," minutes.")' . ,(end_time-stda_time)/60.0 write(*,*) - !CALL EXIT([STATUS]) + !CALL EXIT([STATUS]) endif - + deallocate(apb,ambsqr) !******************************************************************************** open(unit=53,file='amb',form='unformatted',status='old') @@ -890,33 +908,33 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, c big arrays not needed anymore deallocate(pij,qab,pia,ed) - if(.not.printexciton)deallocate(qia) + if(.not.printexciton)deallocate(qia) write(*,*)'diagonalizing ...' write(*,'('' estimated time (min) '',f8.2)') . float(nci)**2*float(nroot)/8.d+8/60. -c if LAPACK does not work +c if LAPACK does not work c allocate(eci(nci),uci(nci,nroot), c . stat=ierr) c call sHQRII(hci,nci,nroot,eci,uci) c faster by a factor of 2-3 - lwork =26*nci + lwork =26*nci liwork=10*nci -c we allocate uci with nci (and not with nroot) as save +c we allocate uci with nci (and not with nroot) as save c choice (other values gave segfaults) nroot=min(nci,int(1.5*nroot)) allocate(eci(nci),uci(nci,nci),work(lwork), . iwork(liwork),isuppz(nci)) if(ierr.ne.0)stop 'allocation failed for TDA matrix diag' - vl=0 + vl=0 vu=thr call ssyevr('V','V','U',nci,hci,nci,vl,vu,il,iu,1.e-6, . nfound,eci,uci,nci,isuppz, . work,lwork,iwork,liwork,info) - + nroot=nfound if(nfound.lt.1) stop 'internal error in diag' @@ -935,7 +953,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ! enddo !!!! - endif + endif write(*,'(i5,'' roots found, lowest/highest eigenvalue : '', .2F8.3,i4)')nroot,eci(1)*27.21139,eci(nroot)*27.21139,info @@ -943,7 +961,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, allocate(umerk(14,nroot)) c contract WF with MO integrals for intensities - p23= 2.0d0/3.0d0 + p23= 2.0d0/3.0d0 aksqrt=sqrt(ak) ! if aniso, store and print the x,y,z-resolved data @@ -953,11 +971,11 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, umrkx=0.0d0 umrky=0.0d0 umrkz=0.0d0 - endif + endif ! if desired, print out data for exciton coupling, i.e., transition moments if(printexciton) then - open(unit=27,file='tda.exc') + open(unit=27,file='tda.exc') write(27,*)ncent,nroot write(27,'(a)',advance='yes')'$coord' ihilf=0 @@ -978,21 +996,21 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do i=1,nroot de=dble(eci(i)) ef=1.0d0/(de+1.0d-8) - xlu=0.0d0 - ylu=0.0d0 - zlu=0.0d0 - xmu=0.0d0 - ymu=0.0d0 - zmu=0.0d0 - xvu=0.0d0 - yvu=0.0d0 + xlu=0.0d0 + ylu=0.0d0 + zlu=0.0d0 + xmu=0.0d0 + ymu=0.0d0 + zmu=0.0d0 + xvu=0.0d0 + yvu=0.0d0 zvu=0.0d0 xms=0.0d0 yms=0.0d0 - zms=0.0d0 + zms=0.0d0 l=0 - umax=-1 - kmem=1 + umax=-1 + kmem=1 q1=0.0 c loop over CSF do j=1,nci @@ -1023,7 +1041,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ymu=ymu+ym(ij)*pp zmu=zmu+zm(ij)*pp if(.not.velcorr.and.printexciton) then - l=(io-1)*nv+(iv-no) + l=(io-1)*nv+(iv-no) q1(1:ncent)=q1(1:ncent)+qia(1:ncent,l)*real(uu) endif enddo @@ -1039,8 +1057,8 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, ymu=ymu*aksqrt zmu=zmu*aksqrt -c polarizability - alp_real(1)=alp_real(1)+xlu*xlu*ef +c polarizability + alp_real(1)=alp_real(1)+xlu*xlu*ef alp_real(2)=alp_real(2)+xlu*ylu*ef alp_real(3)=alp_real(3)+ylu*ylu*ef alp_real(4)=alp_real(4)+xlu*zlu*ef @@ -1059,20 +1077,20 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do k=1,ncent write(27,'(f12.8)') real(aksqrt)*q1(k) ! factor of 2**0.5 arises from spin integration enddo -! write(27,*) +! write(27,*) endif - xp=xlu*xlu+ylu*ylu+zlu*zlu - fl=de*p23*xp - xp=xvu*xvu+yvu*yvu+zvu*zvu - fv=ef*p23*xp + xp=xlu*xlu+ylu*ylu+zlu*zlu + fl=de*p23*xp + xp=xvu*xvu+yvu*yvu+zvu*zvu + fv=ef*p23*xp xp=xlu*xmu+ylu*ymu+zlu*zmu - rl=-235.7220d0*xp -! multiplying with ak (2=singlet excitation -> restricted case) is already included in moments: see book by Harada & Nakanishi + rl=-235.7220d0*xp +! multiplying with ak (2=singlet excitation -> restricted case) is already included in moments: see book by Harada & Nakanishi ! recalculated value for Rot with recent CODATA values and changed it (used to be 235.730) ! the factor is the conversion from e*a0*mu_b in atomic units to 10^40 cgs (m_b is the Bohr magneton and a0 the Bohr radius) - xp=xvu*xmu+yvu*ymu+zvu*zmu - rv=-235.7220d0*xp*ef + xp=xvu*xmu+yvu*ymu+zvu*zmu + rv=-235.7220d0*xp*ef rvp(i)=rv c sum rule sumf=sumf+fl @@ -1148,39 +1166,39 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, endif enddo - + c sort (output) c do 140 ii = 2,nroot -c i = ii - 1 -c k = i -c pp= umerk(1,i) -c do 120 j = ii, nroot -c if (umerk(1,j) .ge. pp) go to 120 -c k = j -c pp=umerk(1,j) -c 120 continue -c if (k .eq. i) go to 140 -c umerk(1,k) = umerk(1,i) -c umerk(1,i) = pp +c i = ii - 1 +c k = i +c pp= umerk(1,i) +c do 120 j = ii, nroot +c if (umerk(1,j) .ge. pp) go to 120 +c k = j +c pp=umerk(1,j) +c 120 continue +c if (k .eq. i) go to 140 +c umerk(1,k) = umerk(1,i) +c umerk(1,i) = pp c do m=2,11 -c hilf=umerk(m,i) -c umerk(m,i)=umerk(m,k) -c umerk(m,k)=hilf +c hilf=umerk(m,i) +c umerk(m,i)=umerk(m,k) +c umerk(m,k)=hilf c enddo -c 140 continue +c 140 continue 11 format(i5,f9.3,f8.1,2f11.4,3x,3(F6.2,'(',i4,'->',i4,')')) 12 format(i5,f6.2,f8.1, 5x,i4,' ->',i4,5x,'gap,J,K:',3f8.3, . 3x,'Kshft:',f8.3) 13 format( 41x,f8.3,13x,f8.3) - + if(.not.rpachk) deallocate(hci) if(.not.velcorr) then write(*,*)'writing spectral data to tda.dat ...' - call print_tdadat(nroot,xmolw,eci,umerk(2:5,:),thr,'tda.dat') + call print_tdadat(nroot,xmolw,eci,umerk(2:5,:),thr,'tda.dat') if(aniso) then write(*,*)'writing anisotropic data ...' @@ -1192,40 +1210,40 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, call print_tdadat(nroot,xmolw,eci,umrkz(1:4,:),thr,'tdaz.dat') deallocate(umrkx,umrky,umrkz) write(*,*)'data given such that f = (f_x + f_y + f_z)/3' - write(*,*)' and R = R_x + R_y + R_z' - endif + write(*,*)' and R = R_x + R_y + R_z' + endif else - ! if velocity correction is on, get improved rotatory strengths - if (printexciton) then + ! if velocity correction is on, get improved rotatory strengths + if (printexciton) then call apbtrafoexc(nci,nroot,uci,eci,xl,yl,zl,xv,yv,zv,xm,ym,zm . ,xmolw,no,nv,coc,ncent,qia,maxconf,iconf,ak,rvp) else call apbtrafo(nci,nroot,uci,eci,xl,yl,zl,xv,yv,zv,xm,ym,zm - . ,xmolw,maxconf,iconf,ak,rvp) + . ,xmolw,maxconf,iconf,ak,rvp) endif do i=1,nroot umerk(5,i)=rvp(i) - enddo + enddo endif - if(printexciton) then + if(printexciton) then deallocate(qia) close(27) endif - -c output + +c output write(*,*) write(*,*) .'excitation energies, transition moments and TDA amplitudes' - if(velcorr) then + if(velcorr) then write(*,*) 'state eV nm fL Rv(corr)' else write(*,*) 'state eV nm fL Rv' endif do i=1,nroot ec=umerk(1,i) - write(*,11) + write(*,11) . i,ec*27.21139, . 1.d+7/(ec*2.19474625d+5),umerk(2,i),umerk(5,i) . ,umerk(6,i),idint(umerk(7,i)),idint(umerk(8,i)) @@ -1241,10 +1259,10 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, call sosor(nroot,xmolw,eci,rvp) deallocate(rvp,q1) - + ***************************** c Lin. Response func. ESA * -***************************** +***************************** if(ESA)then if(triplet) stop 'not available' ! Unrelaxed state-to-state transition dipole moments @@ -1267,9 +1285,9 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, . ,(end_time-stda_time)/60.0 endif write(*,*) - !CALL EXIT([STATUS]) + !CALL EXIT([STATUS]) endif - + ! call hyperpol_sos(nroot,nci,eci,uci,hci,xl,yl,zl,moci, ! .maxconf,iconf,no,nv) !to test the moments ! call tpa_sos(nroot,nci,eci,uci,hci,xl,yl,zl,moci, @@ -1301,13 +1319,13 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, call printvectda(rpachk,nci,nroot,uci,eci) ! TDA vectors endif endif - + if(nto)then if(rpachk)then call print_nto_rpa(uci,hci,ca,moci,nci,nroot,nao,iconf, . maxconf,no,nv) else - + call print_nto(uci,ca,moci,nci,nroot,nao,iconf,maxconf,no,nv) endif endif @@ -1316,44 +1334,44 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, if(rpachk) then write(*,*) 'sTD-DFT done.' - else + else write(*,*) 'sTDA done.' endif !!! old fitting part -!c used for fitting +!c used for fitting inquire(file='.REF',exist=ex) if(.not.ex) return open(unit=27,file='.REF') open(unit=28,file='.OUT') read(27,*) i,k - + if(k.eq.1)then -c take only bright ones +c take only bright ones read(27,*) hilf do j=1,nroot if(umerk(2,j).gt.0.03)then - write(28,'(4f12.6)') + write(28,'(4f12.6)') . (umerk(1,j)*27.211-hilf),hilf,umerk(1,j)*27.211,umerk(2,j) exit endif enddo elseif(k.eq.2)then -c take only very bright ones +c take only very bright ones read(27,*) hilf do j=1,nroot if(umerk(2,j).gt.0.2)then - write(28,'(4f12.6)') + write(28,'(4f12.6)') . (umerk(1,j)*27.211-hilf),hilf,umerk(1,j)*27.211,umerk(2,j) exit endif enddo elseif(k.eq.3)then -c take only very bright CD ones +c take only very bright CD ones read(27,*) hilf do j=1,nroot if(abs(umerk(5,j)).gt.100)then - write(28,'(4f12.6)') + write(28,'(4f12.6)') . (umerk(1,j)*27.211-hilf),hilf,umerk(1,j)*27.211,umerk(5,j) exit endif @@ -1362,7 +1380,7 @@ SUBROUTINE stda(ncent,nmo,nao,xyz,c,eps,occ,iaoat,thr,thrp, do j=1,i read(27,*) hilf if(hilf.gt.0.01) - . write(28,'(4f12.6)') + . write(28,'(4f12.6)') . (umerk(1,j)*27.211-hilf),hilf,umerk(1,j)*27.211,umerk(2,j) enddo endif @@ -1404,12 +1422,12 @@ subroutine makel(nao, s, x) enddo call dsyev ('V','U',nao,x,nao,e,aux,lwork,info) - + c call dHQRII(s,nao,nao,e,vecs) do i=1,nao - - if(bfw)then + + if(bfw)then if(e(i).lt.0)then write(*,*) '!!!!! BIG FAT WARNING !!!!!' write(*,*) '!!!!! BIG FAT WARNING !!!!!' @@ -1442,7 +1460,7 @@ subroutine makel(nao, s, x) enddo !$omp parallel private (i,m) !$omp do - do m=1,nao + do m=1,nao do i=1,nao vecs(i,m)= x(i,m) cc(i,m)=e(m)*x(i,m) @@ -1490,22 +1508,22 @@ subroutine shrink4(n,x,y) end *********************************************************************** -* atomic Löwdin population for a (transition) density with MOs ij +* atomic Löwdin population for a (transition) density with MOs ij *********************************************************************** subroutine lo12pop(i,j,ncent,nao,iaoat,ca,q) - IMPLICIT NONE + IMPLICIT NONE integer i,j,nao,ncent,iaoat(*) real*8 ca(*) real*4 q(*) - integer k,ii,jj,iat - + integer k,ii,jj,iat + q(1:ncent) = 0.0 ii=(i-1)*nao jj=(j-1)*nao - do k=1,nao + do k=1,nao iat=iaoat(k) q(iat)=q(iat)+ca(k+ii)*ca(k+jj) enddo @@ -1522,13 +1540,13 @@ subroutine kshift_to_ediag(de,ek) real*8, intent(inout) :: de ! A(ia,ia) real*8, intent (in) :: ek ! K(ia,ia) real*8 shft,wau,sau,ekia,dmp - ! we apply a rationally damped, K(ia,ia) dependent shift to A(ia,ia) + ! we apply a rationally damped, K(ia,ia) dependent shift to A(ia,ia) ! ! f [ K(ia,ia) ] = shiftmax /( 1 + (K(ia,ia)/width)**steepness ! ! A(ia,ia)=A(ia,ia) + f [ K(ia,ia) ] - ! shiftmax and width are given eV, convert to a.u. + ! shiftmax and width are given eV, convert to a.u. sau=shftmax*0.036749320d0 wau=shftwidth*0.036749320d0 shft = 1.0d0 + (ek/wau)**shftsteep @@ -1559,7 +1577,7 @@ subroutine ptselect(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, logical, allocatable :: incl_conf(:) ak=real(dak) ax=real(dax) - allocate(qj(ncent),qk(ncent),pt(nex),pt2(nex),incl_conf(nexpt), + allocate(qj(ncent),qk(ncent),pt(nex),pt2(nex),incl_conf(nexpt), . stat=ierr) if(ierr.ne.0)stop 'allocation for PT intermediates crashed' incl_conf=.false. @@ -1571,7 +1589,7 @@ subroutine ptselect(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, ! loop over secondary/neglected CSF, this is done omp-parallel !$omp parallel private(i,k,io,iv,iiv,iwrk,qk,de,j,jo,jv,jjv,jwrk,ek,ej, !$omp& qj,pt,amat,pert) reduction (+:pt2) -!$omp do +!$omp do do k=1,nexpt io=kconf(k,1) iv=kconf(k,2) @@ -1579,7 +1597,7 @@ subroutine ptselect(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, iwrk=(io-1)*nv + iiv qk(1:ncent)=pia(1:ncent,iwrk) de=edpt(k) -c loop over primary CSF +c loop over primary CSF do j=1,nex jo=iconf(j,1) jv=iconf(j,2) @@ -1590,16 +1608,16 @@ subroutine ptselect(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, qj(1:ncent)=pij(1:ncent,lin(io,jo)) ej=sdot(ncent,qj,1,qab(1,lin(iiv,jjv)),1) pt(j)=0.0d0 - amat=ak*ek-ej + amat=ak*ek-ej pt(j)=amat**2/(de-ed(j)+1.d-10) enddo pert=sum(pt(1:nex)) -c if sum > threshold include it +c if sum > threshold include it if(pert.gt.thrp)then incl_conf(k)=.true. else do i=1,nex - pt2(i)=pt2(i)+pt(i) + pt2(i)=pt2(i)+pt(i) enddo endif enddo @@ -1609,19 +1627,19 @@ subroutine ptselect(nex,ncent,no,nv,nexpt,mxcnf,iconf,kconf, do i=1,nexpt if(.not.incl_conf(i)) cycle io=kconf(i,1) - iv=kconf(i,2) - de=edpt(i) + iv=kconf(i,2) + de=edpt(i) new=new+1 - iconf(nex+new,1)=io - iconf(nex+new,2)=iv + iconf(nex+new,1)=io + iconf(nex+new,2)=iv ed (nex+new )=de enddo pert=0.0d0 do i=1,nex - pert=pert+pt2(i) - ed(i)=ed(i)-pt2(i) + pert=pert+pt2(i) + ed(i)=ed(i)-pt2(i) enddo - + write(*,'('' average/max PT2 energy lowering (eV):'',2F10.3)') . 27.21139*pert/float(nex),maxval(pt2)*27.21139 @@ -1648,7 +1666,7 @@ subroutine rrpamat(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, real*8, intent(in) :: dak,dax,ed(mxcnf) real*4, intent(out) :: apb(nci*(nci+1)/2),ambsqr(nci*(nci+1)/2) real*4, allocatable :: qk(:),qj(:) - integer i,j,ij,io,iv,jo,jv,ierr,lin,iiv,jjv,iwrk,jwrk + integer i,j,ij,io,iv,jo,jv,ierr,lin,iiv,jjv,iwrk,jwrk real*4 ek,ej,sdot,ak,ax,de allocate(qk(ncent), stat=ierr) if(ierr.ne.0)stop 'allocation for q1 crashed' @@ -1678,17 +1696,17 @@ subroutine rrpamat(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, ek=sdot(ncent,qk,1,qia(1,jwrk),1) ! ek = (ia|jb) apb(ij)=2.0*ak*ek ambsqr(ij)=0.0 - enddo ! j + enddo ! j de=real(ed(i)) ij=lin(i,i) - ek=sdot(ncent,qk,1,qia(1,iwrk),1) + ek=sdot(ncent,qk,1,qia(1,iwrk),1) ambsqr(ij)=sqrt(de-ak*ek) ! diagonal element of (A-B)^0.5 apb(ij)=de+ak*ek ! diagonal element of A+B enddo ! i !$omp end do !$omp end parallel - else - allocate(qj(ncent), stat=ierr) + else + allocate(qj(ncent), stat=ierr) ij=0 ! for now ambsqr=A+B and apb=A-B, since we need to take the sqrt of A-B (but want to save memory) !$omp parallel private(ij,i,j,io,iv,jo,jv,iiv,iwrk,jjv,jwrk,qk,qj,ek,ej) @@ -1729,28 +1747,28 @@ subroutine rrpamat(nci,ncent,no,nv,mxcnf,iconf,dak,dax,ed, ! call prmat4(6,apb,nci,0,'A-B') ! call prmat4(6,ambsqr,nci,0,'A+B') - if(aresp.or.resp) then + if(aresp.or.resp.or.optrota) then write(*,*) ' calculating (A-B)^0.5 not necessary...' - open(unit=53,file='amb',form='unformatted',status='replace') + open(unit=53,file='amb',form='unformatted',status='replace') write(53) apb close(53) apb=ambsqr else - open(unit=52,file='apbmat',form='unformatted',status='replace') + open(unit=52,file='apbmat',form='unformatted',status='replace') write(52) ambsqr - open(unit=53,file='amb',form='unformatted',status='replace') + open(unit=53,file='amb',form='unformatted',status='replace') write(53) apb write(*,*) ' calculating (A-B)^0.5 ...' write(*,'('' estimated time (min) '',f8.2)') . float(nci)**2*float(nci)/4.d+8/60. - call smatpow(nci,apb) ! calculate sqrt(a-b), termed ambsqr + call smatpow(nci,apb) ! calculate sqrt(a-b), termed ambsqr ambsqr=apb rewind(52) read(52) apb close(52,status='delete') - close(53) + close(53) endif deallocate(qj) @@ -1804,10 +1822,10 @@ subroutine rtdamat(nci,ncent,no,nv,mxcnf,iconf,dak,dax hci(j,i)=ak*ek-ej hci(i,j)=hci(j,i) enddo - hci(i,i)=real(ed(i)) + hci(i,i)=real(ed(i)) enddo !$omp end do -!$omp end parallel +!$omp end parallel deallocate(qk,qj) return end subroutine rtdamat @@ -1825,107 +1843,107 @@ subroutine setrep(rep) c values in the paper multiplied by two because (ii:ii)=(IP-EA)=d^2 E/dN^2 but the hardness c definition they use is 1/2d^2 E/dN^2 (in Eh) - rep( 1)= 0.472592880d0 - rep( 2)= 0.922033910d0 - rep( 3)= 0.174528880d0 - rep( 4)= 0.257007330d0 - rep( 5)= 0.339490860d0 - rep( 6)= 0.421954120d0 - rep( 7)= 0.504381930d0 - rep( 8)= 0.586918630d0 - rep( 9)= 0.669313510d0 - rep(10)= 0.751916070d0 - rep(11)= 0.179641050d0 - rep(12)= 0.221572760d0 - rep(13)= 0.263485780d0 - rep(14)= 0.305396450d0 - rep(15)= 0.347340140d0 - rep(16)= 0.389247250d0 - rep(17)= 0.431156700d0 - rep(18)= 0.473082690d0 - rep(19)= 0.171054690d0 - rep(20)= 0.202762440d0 - rep(21)= 0.210073220d0 - rep(22)= 0.217396470d0 - rep(23)= 0.224710390d0 - rep(24)= 0.232015010d0 - rep(25)= 0.239339690d0 - rep(26)= 0.246656380d0 - rep(27)= 0.253982550d0 - rep(28)= 0.261288630d0 - rep(29)= 0.268594760d0 - rep(30)= 0.275925650d0 - rep(31)= 0.307629990d0 - rep(32)= 0.339315800d0 - rep(33)= 0.372359850d0 - rep(34)= 0.402735490d0 - rep(35)= 0.434457760d0 - rep(36)= 0.466117080d0 - rep(37)= 0.155850790d0 - rep(38)= 0.186493240d0 - rep(39)= 0.193562100d0 - rep(40)= 0.200633110d0 - rep(41)= 0.207705220d0 - rep(42)= 0.214772540d0 - rep(43)= 0.221846140d0 - rep(44)= 0.228918720d0 - rep(45)= 0.235986210d0 - rep(46)= 0.243056120d0 - rep(47)= 0.250130180d0 - rep(48)= 0.257199370d0 - rep(49)= 0.287847800d0 - rep(50)= 0.318486730d0 - rep(51)= 0.349124310d0 - rep(52)= 0.379765930d0 - rep(53)= 0.410408080d0 - rep(54)= 0.441057770d0 - rep(55)= 0.050193320d0 - rep(56)= 0.067625700d0 - rep(57)= 0.085044450d0 - rep(58)= 0.102477360d0 - rep(59)= 0.119911050d0 - rep(60)= 0.137327720d0 - rep(61)= 0.154762970d0 - rep(62)= 0.172182650d0 - rep(63)= 0.189612880d0 - rep(64)= 0.207047600d0 - rep(65)= 0.224467520d0 - rep(66)= 0.241896450d0 - rep(67)= 0.259325030d0 - rep(68)= 0.276760940d0 - rep(69)= 0.294182310d0 - rep(70)= 0.311595870d0 - rep(71)= 0.329022740d0 - rep(72)= 0.345922980d0 - rep(73)= 0.363880480d0 - rep(74)= 0.381305860d0 - rep(75)= 0.398774760d0 - rep(76)= 0.416142980d0 - rep(77)= 0.433645100d0 - rep(78)= 0.451040140d0 - rep(79)= 0.468489860d0 - rep(80)= 0.485845500d0 - rep(81)= 0.125267300d0 - rep(82)= 0.142686770d0 - rep(83)= 0.160116150d0 - rep(84)= 0.177558890d0 - rep(85)= 0.194975570d0 - rep(86)= 0.212407780d0 - rep(87)= 0.072635250d0 - rep(88)= 0.094221580d0 - rep(89)= 0.099202950d0 - rep(90)= 0.104186210d0 - rep(91)= 0.142356330d0 - rep(92)= 0.163942940d0 - rep(93)= 0.185519410d0 - rep(94)= 0.223701390d0 + rep( 1)= 0.472592880d0 + rep( 2)= 0.922033910d0 + rep( 3)= 0.174528880d0 + rep( 4)= 0.257007330d0 + rep( 5)= 0.339490860d0 + rep( 6)= 0.421954120d0 + rep( 7)= 0.504381930d0 + rep( 8)= 0.586918630d0 + rep( 9)= 0.669313510d0 + rep(10)= 0.751916070d0 + rep(11)= 0.179641050d0 + rep(12)= 0.221572760d0 + rep(13)= 0.263485780d0 + rep(14)= 0.305396450d0 + rep(15)= 0.347340140d0 + rep(16)= 0.389247250d0 + rep(17)= 0.431156700d0 + rep(18)= 0.473082690d0 + rep(19)= 0.171054690d0 + rep(20)= 0.202762440d0 + rep(21)= 0.210073220d0 + rep(22)= 0.217396470d0 + rep(23)= 0.224710390d0 + rep(24)= 0.232015010d0 + rep(25)= 0.239339690d0 + rep(26)= 0.246656380d0 + rep(27)= 0.253982550d0 + rep(28)= 0.261288630d0 + rep(29)= 0.268594760d0 + rep(30)= 0.275925650d0 + rep(31)= 0.307629990d0 + rep(32)= 0.339315800d0 + rep(33)= 0.372359850d0 + rep(34)= 0.402735490d0 + rep(35)= 0.434457760d0 + rep(36)= 0.466117080d0 + rep(37)= 0.155850790d0 + rep(38)= 0.186493240d0 + rep(39)= 0.193562100d0 + rep(40)= 0.200633110d0 + rep(41)= 0.207705220d0 + rep(42)= 0.214772540d0 + rep(43)= 0.221846140d0 + rep(44)= 0.228918720d0 + rep(45)= 0.235986210d0 + rep(46)= 0.243056120d0 + rep(47)= 0.250130180d0 + rep(48)= 0.257199370d0 + rep(49)= 0.287847800d0 + rep(50)= 0.318486730d0 + rep(51)= 0.349124310d0 + rep(52)= 0.379765930d0 + rep(53)= 0.410408080d0 + rep(54)= 0.441057770d0 + rep(55)= 0.050193320d0 + rep(56)= 0.067625700d0 + rep(57)= 0.085044450d0 + rep(58)= 0.102477360d0 + rep(59)= 0.119911050d0 + rep(60)= 0.137327720d0 + rep(61)= 0.154762970d0 + rep(62)= 0.172182650d0 + rep(63)= 0.189612880d0 + rep(64)= 0.207047600d0 + rep(65)= 0.224467520d0 + rep(66)= 0.241896450d0 + rep(67)= 0.259325030d0 + rep(68)= 0.276760940d0 + rep(69)= 0.294182310d0 + rep(70)= 0.311595870d0 + rep(71)= 0.329022740d0 + rep(72)= 0.345922980d0 + rep(73)= 0.363880480d0 + rep(74)= 0.381305860d0 + rep(75)= 0.398774760d0 + rep(76)= 0.416142980d0 + rep(77)= 0.433645100d0 + rep(78)= 0.451040140d0 + rep(79)= 0.468489860d0 + rep(80)= 0.485845500d0 + rep(81)= 0.125267300d0 + rep(82)= 0.142686770d0 + rep(83)= 0.160116150d0 + rep(84)= 0.177558890d0 + rep(85)= 0.194975570d0 + rep(86)= 0.212407780d0 + rep(87)= 0.072635250d0 + rep(88)= 0.094221580d0 + rep(89)= 0.099202950d0 + rep(90)= 0.104186210d0 + rep(91)= 0.142356330d0 + rep(92)= 0.163942940d0 + rep(93)= 0.185519410d0 + rep(94)= 0.223701390d0 return end **************************************************************************** -* Given orbital indeces i1 and i2, lin() returns index in the linear array * +* Given orbital indeces i1 and i2, lin() returns index in the linear array * **************************************************************************** integer*4 function lin(i1,i2) @@ -1936,7 +1954,7 @@ integer*4 function lin(i1,i2) lin=idum2+idum1*(idum1-1)/2 return end - + integer*8 function lin8(i1,i2) integer i1,i2 integer*8 idum1,idum2 @@ -1955,24 +1973,24 @@ subroutine cofc(nat,xyz,coc) integer iat,iatyp real*8 ctot ! real*8 atwt(111) -! data atwt / 1.00794, 4.002602, 6.941, 9.012182, 10.811, 12.0107, -! . 14.0067, 15.9994, 18.9984032, 20.1797, 22.98976928, -! . 24.305, 26.9815386, 28.0855, 30.973762, 32.065, -! . 35.453, 39.948, 39.0983, 40.078, 44.955912, 47.867, -! . 50.9415, 51.9961, 54.938045, 55.845, 58.933195, -! . 58.6934, 63.546, 65.38, 69.723, 72.64, 74.9216, -! . 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585, -! . 91.224, 92.90638, 95.96, 98.0, 101.07, 102.90550, -! . 106.42, 107.8682, 112.411, 114.818, 118.71, 121.76, -! . 127.6, 126.90447, 131.293, 132.9054519, 137.327, +! data atwt / 1.00794, 4.002602, 6.941, 9.012182, 10.811, 12.0107, +! . 14.0067, 15.9994, 18.9984032, 20.1797, 22.98976928, +! . 24.305, 26.9815386, 28.0855, 30.973762, 32.065, +! . 35.453, 39.948, 39.0983, 40.078, 44.955912, 47.867, +! . 50.9415, 51.9961, 54.938045, 55.845, 58.933195, +! . 58.6934, 63.546, 65.38, 69.723, 72.64, 74.9216, +! . 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585, +! . 91.224, 92.90638, 95.96, 98.0, 101.07, 102.90550, +! . 106.42, 107.8682, 112.411, 114.818, 118.71, 121.76, +! . 127.6, 126.90447, 131.293, 132.9054519, 137.327, ! . 138.90547, 140.116, 140.90765, 144.242, 145.0, 150.36, -! . 151.964, 157.25, 158.92535, 162.5, 164.93032, -! . 167.259, 168.93421, 173.054, 174.9668, 178.49, +! . 151.964, 157.25, 158.92535, 162.5, 164.93032, +! . 167.259, 168.93421, 173.054, 174.9668, 178.49, ! . 180.94788, 183.84, 186.207, 190.23, 192.217, 195.084, ! . 196.966569, 200.59, 204.3833, 207.2, 208.9804, 209.0, -! . 210.0, 222.0, 223.0, 226.0, 227.0, 232.03806, -! . 231.03588, 238.02891, 237.0, 244.0, 243.0, 247.0, -! . 247.0, 251.0, 252.0, 257.0, 258.0, 259.0, 262.0, +! . 210.0, 222.0, 223.0, 226.0, 227.0, 232.03806, +! . 231.03588, 238.02891, 237.0, 244.0, 243.0, 247.0, +! . 247.0, 251.0, 252.0, 257.0, 258.0, 259.0, 262.0, ! . 267.0, 268.0, 271.0, 272.0, 270.0, 276.0, 281.0, 280.0 / iatyp=0 @@ -1989,19 +2007,19 @@ subroutine cofc(nat,xyz,coc) end subroutine cofc - subroutine print_tdadat(nroot,xmolw,energy,trstr,thr,fname) + subroutine print_tdadat(nroot,xmolw,energy,trstr,thr,fname) implicit none integer, intent( in ) :: nroot real*4, intent( in ) :: energy(*) real*8, intent( in ) :: trstr(4,nroot),thr,xmolw character*80, intent( in ) :: fname integer i - + open(unit=26,file=trim(fname)) write(26,*)'NM' - write(26,*)'UV' + write(26,*)'UV' write(26,*)'MMASS' - write(26,*)xmolw + write(26,*)xmolw write(26,*)'LFAKTOR' write(26,*)' 0.5' write(26,*)'RFAKTOR' @@ -2017,6 +2035,6 @@ subroutine print_tdadat(nroot,xmolw,energy,trstr,thr,fname) . i,energy(i)*27.21139,trstr(1,i),trstr(2,i),trstr(3,i), . trstr(4,i) enddo - close (26) + close (26) - end subroutine print_tdadat + end subroutine print_tdadat diff --git a/stdacommon.f90 b/stdacommon.f90 index 7a38a9f..f8fee42 100644 --- a/stdacommon.f90 +++ b/stdacommon.f90 @@ -15,9 +15,9 @@ ! You should have received a copy of the GNU Lesser General Public License ! along with stda. If not, see . ! -! wavefunction and basis common block +! wavefunction and basis common block module stdacommon - + real*8, allocatable :: co(:,:),exip(:),occ(:),eps(:) real*8, allocatable :: cxip(:),eta(:,:) integer, allocatable :: ipat(:),ipty(:),ipao(:),iaoat(:) @@ -25,7 +25,7 @@ module stdacommon end module stdacommon -! xtb common block +! xtb common block module kshiftcommon real*8 shftmax,shftwidth,shftsteep,shftmax_somo @@ -39,9 +39,10 @@ end module kshiftcommon module commonlogicals logical triplet,rpachk,eigvec,screen,dokshift,printexciton,velcorr logical aniso - logical resp,TPA,aresp,ESA,smp2,bfw,spinflip,rw,pt_off,nto + logical resp,TPA,aresp,ESA,smp2,bfw,spinflip,rw,pt_off,nto,sf_s2 + logical optrota,velo_OR end module commonlogicals - + ! some variables for the response functions module commonresp integer :: num_freq