program cpropag implicit none integer*4 m1,m2,N,nat,natu,nu,n3,nu3, IERR,iq,iwr,it,ix, 1I,itt,umin,umax,np,is,ia,nd real*8 c1,sa,v(3),am(3,3),t1(3),pi,po,ds,rs,at(3,3), 1wmin,wmax,delq,rq,qabs,qvcd,qvmi,qvma,ecm,qmax,wsmin,wsmax,fwhm, 1tau,DSX(3),QSX(3),MSX(3),amn(3,3),ram(13),ENM,temp,vmin,vmax,gau, 1vau(4),ux(4),uy(4),epsr real*8,allocatable::r(:),m(:),f(:,:),Dqc(:,:),Sqr(:,:),Eq(:), 1Dqs(:,:),Sqi(:,:),fv1(:),fv2(:),fm1(:,:),mu(:),P(:,:,:),A(:,:,:), 1sd(:,:),umr(:,:),umi(:,:),qq(:,:),prc(:,:),pic(:,:),arc(:,:), 1aic(:,:),qic(:,:,:),qrc(:,:,:),ALPHA(:,:,:),G(:,:,:),AA(:,:,:,:), 1alrc(:,:,:),alic(:,:,:),grc(:,:,:),gic(:,:,:),aarc(:,:,:,:), 1aaic(:,:,:,:),ru(:),Vr(:,:),Vi(:,:),Cr(:,:),Ci(:,:),Wq(:) integer*4,allocatable::z(:),uu(:,:),iafirst(:),coa(:,:),zu(:) logical LAPT,lsmat,lprn,ltab,lrot,lqpt,ldip,lg,lx,lttt,cpol, 1tdc character*8 QF(8),SF(16) data QF/'DOGX.TAB','DOGY.TAB','DOGZ.TAB',' DOG.TAB', 1 'ROAX.TAB','ROAY.TAB','ROAZ.TAB',' ROA.TAB'/ data SF/' d.prn',' dx.prn',' dy.prn',' dz.prn', 1 ' r.prn',' rx.prn',' ry.prn',' rz.prn', 1 ' ram.prn','ramx.prn','ramy.prn','ramz.prn', 1 ' roa.prn','roax.prn','roay.prn','roaz.prn'/ c1=1302.828d0 pi=4.0d0*atan(1.0d0) po=dsqrt(2.0d0*pi) write(6,601) 601 format(/, 1' Propagate periodic crystal-like 1D structures',/, 2/, 2' Input: FILE.X geometry',/, 1' FILE.F force fieldC',/, 1' FILE.TEN dipole derivatives',/, 1' FILE.QEN quadrupole derivatives',/, 1' FILE.TTT polarizability derivatives',/, 1' FTRY.INP atomic masses',/, 1' CPROPAG.OPT program options',/, 1' ',/, 1' Output: DOG.TAB DOGX.TAB DOGY.TAB DOGZ.TAB',/, 1' - isotropic and oriented IR/VCD intenisites',/, 1' ROA.TAB ROAX.TAB ROAY.TAB ROAZ.TAB',/, 1' - isotropic and oriented RAM/ROA intenisites',/, 1' d.prn dx.prn dy.prn dz,prn',/, 1' - isotropic and oriented IR spectra',/, 1' r.prn rx.prn ry.prn rz,prn',/, 1' - isotropic and oriented VCD spectra',/, 1' ram.prn ramx.prn ramy.prn ramz,prn',/, 1' - isotropic and oriented Raman spectra',/, 1' roa.prn roax.prn roay.prn roaz,prn',/, 1' - isotropic and oriented ROA spectra',/, 1' BIGC.X geometry',/, 1' DIP.TXT transition dipoles',/, 1' Q.PRNi mode contributions',/) allocate(z(1),r(3*1)) CALL READX(0,nat,z,R) n3=3*nat deallocate(z,r) allocate(z(nat),r(n3),uu(nat,nat),m(nat),f(n3,n3),P(nat,3,3), 1A(nat,3,3),qq(3*nat,6),ALPHA(n3,3,3),G(n3,3,3),AA(n3,3,3,3)) CALL READX(1,nat,z,R) CALL READM(nat,m) call readopt(uu,nat,z,nu,natu,m1,m2,N,iwr,LAPT,qmax,lsmat,wmin, 1wmax,delq,umax,umin,qvmi,qvma,ltab,lprn,np,wsmin,wsmax,fwhm,lrot, 1lqpt,ldip,lg,lx,lttt,ENM,temp,cpol,nd,gau,vmin,vmax,vau,ux,uy, 1epsr,is,tdc) c determine first atom in each unit: allocate(iafirst(nu)) call mkfirst(nat,nu,iafirst,uu,natu) nu3=3*natu CALL READFF(n3,f) CALL MASSFF(n3,f,m) CALL READTEN(nat,P,A,LAPT,R) call reqd(lqpt,qq,nat,R,P) if(lttt)then CALL READTTT(ALPHA,AA,G,nat) if(iwr.eq.2)write(6,*)'Common:' if(iwr.eq.2)call writetttd(nat,ALPHA,G,AA) CALL TTTT(n3,ALPHA,AA,G,NAT,1,R) if(iwr.eq.2)write(6,*)'Local:' if(iwr.eq.2)call writetttd(nat,ALPHA,G,AA) endif if(lrot)then c put tensors to common system (only first unit will be used): CALL tenback(nat,P,A,R) call tq(qq,nat,R,P,1) if(lttt)CALL TTTT(n3,ALPHA,AA,G,NAT,2,R) c produce Cartesian system rotation matrix, shift vector and twist: c at - first index old (input) system call crmi(r,nat,natu,uu,m1,m2,v,at,T1,tau) allocate(umr(3,3),umi(3,3),prc(nu3,3),pic(nu3,3),arc(nu3,3), 1 aic(nu3,3),qrc(nu3,3,3),qic(nu3,3,3)) c rotate FF and tensors to the propagation Cartesian system, c propagation axis = z: call rot1(f,at,n3) call rot2(P,at,n3) call rot2(A,at,n3) call rot3(qq,at,n3) if(lttt)then call rot4(ALPHA,at,n3) call rot4(G,at,n3) call rot5(AA,at,n3) endif c make the matrix umr umi into complex coordinates: call setum(umr,umi) c rotate tensors into complex coordinates, c first cell atoms only, transform dipole indices only: call rott(P ,prc,pic,nat,umr,umi,natu,uu,m1) call rott(A ,arc,aic,nat,umr,umi,natu,uu,m1) call rotQ(qq,qrc,qic,nat,umr,umi,natu,uu,m1) allocate(alrc(nu3,3,3),alic(nu3,3,3),grc(nu3,3,3),gic(nu3,3,3), 1 aarc(nu3,3,3,3),aaic(nu3,3,3,3)) if(lttt)then call rottt(ALPHA,alrc,alic,nat,umr,umi,natu,uu,m1) call rottt(G ,grc, gic,nat,umr,umi,natu,uu,m1) call rott3(AA,aarc,aaic,nat,umr,umi,natu,uu,m1) endif c rotation matrix in new coordinates: call srm(amn,tau) c rotation matrix in old coordinates: call rmnc(am,amn,at) else c produce Cartesian rotation matrix and shift vector: call crms(r,nat,natu,uu,m1,m2,v,am,T1) endif allocate(Dqc(nu3,nu3),Sqr(nu3,nu3),Eq(nu3),Dqs(nu3,nu3), 1Sqi(nu3,nu3),fv1(nu3),fv2(nu3),fm1(2,nu3),mu(natu)) c control geometry: if(lx)call gengeo(N,nat,nu3,natu,v,am,r,m1,uu,z,T1) call massa(nat,natu,uu,m,mu,m1) if(iwr.gt.0)write(6,*)' Masses assigned' if(ldip)OPEN(29,FILE='DIP.TXT') if(ltab)then OPEN(18,FILE='DOG.TAB') WRITE (18,501) 501 FORMAT (/,' FREQ. ', 1 ' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]', 2 /,'-----------------------------------------------------------') do 59 i=1,3 open(18+i,file=QF(i)) 59 write(18+i,501) if(lttt)then do 62 i=5,8 call setis(QF(I),is) open(18+i,file=QF(i)(is:len(QF(i)))) 62 write(18+i,501) endif endif if(lprn)then allocate(sd(16,np)) sd=0.0d0 c spectra: 1 absorbtion c 2 abs x c 3 abs y c 4 abs z c 5 vcd c 6 vcd x c 7 vcd y c 8 vcd z c 9 ram c 10 ramx c 11 ramy c 12 ramz c 13 roa c 14 roax c 15 roay c 16 roaz endif if(lsmat)open(46,file='F.INP.SCR') close(46,status='delete') if(lsmat)open(46,file='F.INP.SCR') open(47,file='Q.PRN') ierr=0 it=0 itt=0 allocate(zu(natu),ru(nu3)) c put rel coord of first unit in ru call pa1(m1,nat,natu,z,uu,r,T1,zu,ru) if(cpol)then c correction to group polarizabilities c recognize amides in the first unit call raa(ia,ru,zu,natu) write(6,608)ia 608 format(i6,' amides') allocate(coa(ia,6)) call rda(ia,coa) allocate(Vr(nu3,nu3),Vi(nu3,nu3),Cr(nu3,nu3),Ci(nu3,nu3),Wq(nu3)) endif c rq=-delq write(6,603)0.0d0,qmax,delq 603 format(' Qmin:',f10.3,' Qmax:',f10.3,' DelQ:',f10.3) c loop over wave vectors: c qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq 1000 rq=rq+delq if(rq.gt.qmax)goto 66 write(6,600)rq 600 format(' q =',f10.3) qabs=0.0d0 qvcd=0.0d0 call maked(rq,nat,natu,n3,nu3,nu,f,m1,Dqc,Dqs,N,uu,iwr,iafirst) call lookatme(DQc,DQs,nu3,0) if(tdc)call addtdc(rq,nat,natu,n3,nu3,nu,f,m1,Dqc,Dqs,N,uu, 1ru,am,v,P,iafirst,epsr,nd) c here eigenvalues are force constants, w^2: call ch(nu3,Dqc,Dqs,Eq,1,Sqr,Sqi,fv1,fv2,fm1,ierr) if(ierr.ne.0)call report('diagonalization not successfull') c now sum(Sqr(k,i)*Sqr(k,j)+Sqi(k,i)*Sqi(k,j)) = delta(i,j) c make mass-weighted S-matrix: call sm(Sqr,Sqi,mu,nu3) do 11 iq=1,nu3 sa=1.0d0 if(Eq(iq).lt.0.0d0)sa=-sa 11 Eq(iq)=sa*sqrt(dabs(Eq(iq)))*c1 if(iwr.gt.0)write(6,602)Eq 602 format(6f12.2) call lookate(Eq,nu3,0) c if(cpol)then c correction to group polarizabilities c set dynamic matrix in normal modes call ddm(N,nat,nu3,nd,m1,ia,coa,Sqr,Sqi,Vr,Vi,Eq,ru,P,A, 1 gau,vmin,vmax,rq,uu,am,epsr,vau,v,ux,uy,is) call lookatme(Vr,Vi,nu3,1) write(6,605)nu3,nu3 605 format(' Re-diagonalizing',i7,' x',i6) c here eigenvalues are energies in au, w: call ch(nu3,Vr,Vi,Wq,1,Cr,Ci,fv1,fv2,fm1,ierr) if(ierr.ne.0)call report('diagonalization not successfull') call lookate(Wq,nu3,1) c rewrite old energies, get new effective modes: call tcc(nu3,Cr,Ci,Sqr,Sqi,Eq,Wq) endif c lopp over all modes of this wave vector: c iqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiq do 1 iq=1,nu3 ecm=Eq(iq) if(ecm.gt.1.0d0)then it=it+1 if(lrot)then call getdipr(iq,nu3,prc,pic,arc,aic,qrc,qic,Sqr,Sqi,N,ecm,rq, 1 tau,ds,rs,DSX,MSX,QSX,v,umr,umi,at,ldip,lttt,alrc,alic,grc,gic, 1 aarc,aaic,ENM,ram,iwr) else call getdips(N,nat,natu,nu3,uu,m1,r,iafirst,rq,ecm,lsmat,wmin, 1 wmax,qvmi,qvma,itt,it,umin,umax,am,v,DSX,QSX,MSX,ds,rs, 1 P,qq,A,Sqr,Sqi,qmax,delq,iq,ldip,ALPHA,G,AA,lttt,ram,ENM,iwr) endif if(ltab)then WRITE(18,523)it,ecm,ds,rs,rs,po*Eq(iq)*ds do 60 ix=1,3 60 WRITE(18+ix,523)it,ecm,DSX(ix),QSX(ix)+MSX(ix) 523 FORMAT(I5,f16.4,4G12.4) if(lttt)then do 67 ix=5,7 67 WRITE(18+ix,3000)it,ecm,ram(1),ram(2),ram(3+ix),ram(4),ram(4), 1 ram(6+ix),ram(4),ram(6+ix),ram(6),ram(7) WRITE(18+8,3000)it,ecm,ram(1),ram(2),ram(3),ram(4),ram(4), 1 ram(5),ram(4),ram(5),ram(6),ram(7) 3000 FORMAT(I5,f9.2,6g12.3,f6.3,4g12.4) endif endif if(lprn)then call sspread(ecm,ds,DSX(1),DSX(2),DSX(3),rs,QSX(1)+MSX(1), 1 QSX(2)+MSX(2),QSX(3)+MSX(3),sd,np,wsmin,wsmax,fwhm,lg) if(lttt)call spreadr(ecm,ram,sd,np,wsmin,wsmax,fwhm,lg,temp) endif qabs=qabs+ds qvcd=qvcd+dabs(rs) endif 1 continue c iqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiqiq write(47,4747)rq,qabs,qvcd 4747 format(f12.6,2g12.4) if(rq.lt.qmax)goto 1000 c qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq 66 if(ltab)then do 61 ix=0,3 WRITE(18+ix,1819) 1819 FORMAT('-------------------------------------------------') 61 close(18+ix) write(6,*)'DOG.TAB DOGX.TAB DOGY.TAB DOGZ.TAB written' if(lttt)then do 63 ix=5,8 WRITE(18+ix,1819) 63 close(18+ix) write(6,*)'ROA.TAB ROAX.TAB ROAY.TAB ROAZ.TAB written' endif endif close(47) if(ldip)close(29) if(lprn)then do 64 ix=1,8 call setis(SF(ix),is) 64 call wspectra(SF(ix)(is:len(SF(ix))),wsmin,wsmax,sd,np,ix) if(lttt)then do 65 ix=9,16 call setis(SF(ix),is) 65 call wspectra(SF(ix)(is:len(SF(ix))),wsmin,wsmax,sd,np,ix) endif endif if(lsmat)call smatout(itt,umin,umax,natu,r,N,uu,nat,m1,am,v, 1nu3,z) end subroutine ddm(N,nat,nu3,nd,m1,ia,coa,Sr1,Si1,Dr,Di,e,ru,APT,AAT, 1gau,wmin,wmax,q,uu,A,epsr,vau,v,ux,uy,is) c construct potential matrix between normal modes c I J c D = <0 0 .1 0 ..0 | S T ak T S | 0 0 1 ...0> c N .. number of cells c nu3 .. dimension = 3 x Nat in the cell c ia .. number of amides in one unit c coa .. list of the amide atoms c APT,AAT .. atomic polar and axial tensors c nat .. number of atoms in FILE.X c nd .. units -nd .... nd c ru .. relative coords of atoms of the ref unit c e .. old energies (normal mode frequencies) c q .. wave vector c uu .. cell atom index c v .. the translation vector between cells c vau .. amide vibrational frequencies c ux,uy .. defining projection sof amide I vibrational dipoles c is ... indicator of deuteration implicit none integer*4 N,ia,coa(ia,6),i,j,ii,i6,ai,ix,jx,is,kk, 1ib,uu(nat,nat),m1,nd,nat,nu3,xx,Ie,Ke,ik real*8 e(*),CM,x(3),y(3),z(3),gau,wmin,wmax,A(3,3),con, 1co(3),v(3),rk(3),tik(3,3),q,pi,vau(4), 1Sr1(nu3,nu3),Si1(nu3,nu3),APT(nat,3,3),AAT(nat,3,3), 1Dr(nu3,nu3),Di(nu3,nu3),ri(3),f1,mJ,sJ,cJ,ru(*), 1epsr,ux(4),uy(4),wij complex*16 tt,phase,al(3,3) complex*16,allocatable:: ui(:,:) real*8,allocatable::p(:,:),u(:,:,:),ws(:),Sr(:,:),Si(:,:) complex*16,allocatable::tau(:,:),t0(:,:) c write(6,*)'seting dynamic matrix' allocate(p(2*3*nat,3),ws(nu3)) do 101 i=1,nat do 101 ix=1,3 ii=ix+3*(i-1) do 101 jx=1,3 p( ii,jx)=APT(ia,ix,jx) 101 p(3*nat+ii,jx)=AAT(ia,ix,jx) pi=4.0d0*datan(1.0d0) Dr=0.0d0 Di=0.0d0 CM=219474.630d0 c diagonal elements energy in atomic units: do 1 i=1,nu3 ws(i)=e(i)/CM 1 Dr(i,i)=ws(i) f1=2.0d0*pi*q/dble(N) c S-matrix in atomic units allocate(Sr(nu3,nu3),Si(nu3,nu3)) c con=sqrt(AMU) con=dsqrt(1822.888486209d0) do 102 i=1,nu3 do 102 ix=1,nu3 Sr(i,ix)=Sr1(i,ix)/con 102 Si(i,ix)=Si1(i,ix)/con allocate(tau(nu3,(2*nd+1)*ia*3), 1u((2*nd+1)*ia,4,3),ui(ia,3),t0(nu3,(2*nd+1)*ia*3)) tau=0.0d0 t0=0.0d0 c loop over normal modes: do 8 I=1,nu3 if(ws(i)*CM.le.wmax.and.ws(i)*CM.ge.wmin)then c u = S x P do 111 ii=1,ia do 111 ix=1,3 ui(ii,ix)=dcmplx(0.0d0,0.0d0) c loop over amide atoms in cell 0 do 112 i6=1,6 ai=coa(ii,i6) if(ai.ne.0)then do 113 jx=1,3 c index in FILE.X ib=jx+3*(uu(m1,ai)-1) 113 ui(ii,ix)=ui(ii,ix)+dcmplx(Sr(jx+3*(ai-1),i)*p(ib,ix), 1 Si(jx+3*(ai-1),i)*p(ib,ix)) endif 112 continue 111 ui(ii,ix)=ui(ii,ix)/dcmplx(dsqrt(2.0d0*ws(I))) c loop over amides ik=0 do 9 kk=1,ia c loop over elementary cells do 9 Ke=-nd,nd ik=ik+1 c position of amide kk in cell "0": call asx(rk,ru,coa(kk,1)) c rotate to cell Ke: call rtcj(rk,A,Ke) c position in cell Ke (relative to T1) rk(1)=rk(1)+v(1)*dble(Ke) rk(2)=rk(2)+v(2)*dble(Ke) rk(3)=rk(3)+v(3)*dble(Ke) do 9 ii=1,ia c loop over elementary cells do 9 Ie=-nd,nd if(ii.ne.kk.or.Ie.ne.Ke)then c phase relative to m1: mJ=f1*dble(Ie-m1) sJ=sin(mJ) cJ=cos(mJ) phase=dcmplx(cJ,sJ) call asx(ri,ru,coa(ii,1)) call rtcj(ri,A,Ie) ri(1)=ri(1)+v(1)*dble(Ie) ri(2)=ri(2)+v(2)*dble(Ie) ri(3)=ri(3)+v(3)*dble(Ie) call dtik(tik,ri,rk,epsr) c loop over amide atoms in cell 0 do 902 ix=1,3 xx=ix+3*(ik-1) do 902 jx=1,3 if(Ie.eq.0)t0(I,xx)=t0(I,xx)+ui(ii,jx)*dcmplx(tik(jx,ix)) 902 tau(I,xx)=tau(I,xx)+ui(ii,jx)*dcmplx(tik(jx,ix))*phase endif 9 continue endif 8 continue c set vibrational transition dipole moments for each amide ik=0 do 4 kk=1,ia do 4 Ke=-nd,nd ik=ik+1 c setup local axes x ~ CN, z ~ CN x CO, y ~ z x x call asx(rk,ru,coa(kk,1)) do 3 ix=1,3 x( ix)=ru(ix+3*(coa(kk,4)-1))-rk(ix) 3 co(ix)=ru(ix+3*(coa(kk,3)-1))-rk(ix) call nv(x) call vp(z,x,co) call nv(z) call vp(y,z, x) c rotate to cell Ke: call rtcj(x,A,Ke) call rtcj(y,A,Ke) call rtcj(z,A,Ke) c vibrational dipoles: do 4 i=1,4 do 4 ix=1,3 4 u(ik,i,ix)=ux(I)*x(ix)+uy(i)*y(ix) c loop over normal modes I do 2 I=1,nu3 if(ws(I)*CM.le.wmax.and.ws(I)*CM.ge.wmin)then c loop over normal modes J do 201 J=I+1,nu3 if(ws(J)*CM.le.wmax.and.ws(J)*CM.ge.wmin)then c effective frequency wij=dsqrt(ws(i)*ws(j)) tt=dcmplx(0.0d0) ik=0 do 7 kk=1,ia do 7 Ke=-nd,nd ik=ik+1 c set effective polarizability for this frequency: call seta(wij,ia,is,ik,nd,u,vau,al,gau) do 7 ix=1,3 xx=ix+3*(ik-1) do 7 jx=1,3 c t0* x a x t x exp(i 2 pi q J /N) 7 tt=tt+dconjg(t0(I,xx))*al(ix,jx)*tau(J,jx+3*(ik-1)) Di(I,J)=Di(I,J)-2.0d0*dimag(tt) Dr(I,J)=Dr(I,J)-2.0d0*real(tt) Di(J,I)=-Di(I,J) Dr(J,I)= Dr(I,J) endif 201 continue endif 2 continue return end subroutine rtc2(r,A,J) implicit none real*8 r(3,3),f(3,3),A(3,3) integer*4 jj,J,g c rotate rk in 0 to cell Je: do 901 jj=1,abs(J) if(J.gt.0)then do 1 g=1,3 f(1,g)=A(1,1)*r(1,g)+A(1,2)*r(2,g)+A(1,3)*r(3,g) f(2,g)=A(2,1)*r(1,g)+A(2,2)*r(2,g)+A(2,3)*r(3,g) 1 f(1,g)=A(3,1)*r(1,g)+A(3,2)*r(2,g)+A(3,3)*r(3,g) do 2 g=1,3 r(g,1)=A(1,1)*f(g,1)+A(1,2)*f(g,2)+A(1,3)*f(g,3) r(g,2)=A(2,1)*f(g,1)+A(2,2)*f(g,2)+A(2,3)*f(g,3) 2 r(g,1)=A(3,1)*f(g,1)+A(3,2)*f(g,2)+A(3,3)*f(g,3) else do 3 g=1,3 f(1,g)=A(1,1)*r(1,g)+A(2,1)*r(2,g)+A(3,1)*r(3,g) f(2,g)=A(1,2)*r(1,g)+A(2,2)*r(2,g)+A(3,2)*r(3,g) 3 f(1,g)=A(1,3)*r(1,g)+A(2,3)*r(2,g)+A(3,3)*r(3,g) do 4 g=1,3 r(g,1)=A(1,1)*f(g,1)+A(2,1)*f(g,2)+A(3,1)*f(g,3) r(g,2)=A(1,2)*f(g,1)+A(2,2)*f(g,2)+A(3,2)*f(g,3) 4 r(g,1)=A(1,3)*f(g,1)+A(2,3)*f(g,2)+A(3,3)*f(g,3) endif 901 continue return end subroutine rtcj(rk,A,J) implicit none real*8 rk(3),rf(3),A(3,3) integer*4 jj,J c rotate rk in 0 to cell Je: do 901 jj=1,abs(J) if(J.gt.0)then rf(1)=A(1,1)*rk(1)+A(1,2)*rk(2)+A(1,3)*rk(3) rf(2)=A(2,1)*rk(1)+A(2,2)*rk(2)+A(2,3)*rk(3) rf(1)=A(3,1)*rk(1)+A(3,2)*rk(2)+A(3,3)*rk(3) else rf(1)=A(1,1)*rk(1)+A(2,1)*rk(2)+A(3,1)*rk(3) rf(2)=A(1,2)*rk(1)+A(2,2)*rk(2)+A(3,2)*rk(3) rf(1)=A(1,3)*rk(1)+A(2,3)*rk(2)+A(3,3)*rk(3) endif rk=rf 901 continue return end subroutine dtik(tk,ri,rk,epsr) c distances supplied in Angstroms implicit none integer*4 ix,jx real*8 tk(3,3),ri(3),rk(3),dik,rik(3),epsr,bohr bohr=0.529177d0 do 1 ix=1,3 1 rik(ix)=(ri(ix)-rk(ix))/bohr dik=dsqrt(rik(1)**2+rik(2)**2+rik(3)**2) do 2 ix=1,3 do 2 jx=1,3 2 tk(ix,jx)=3.0d0*rik(ix)*rik(jx) do 3 ix=1,3 3 tk(ix,ix)=tk(ix,ix)-dik**2 do 4 ix=1,3 do 4 jx=1,3 4 tk(ix,jx)=tk(ix,jx)/dik**5/epsr return end subroutine asx(ri,r,ai) integer*4 ai real*8 ri(3),r(*) ri(1)=r(1+3*(ai-1)) ri(2)=r(2+3*(ai-1)) ri(3)=r(3+3*(ai-1)) return end subroutine seta(w,ia,is,ik,nd,u,v,a,gau) implicit none integer*4 is,ix,jx,iv,ik,ia,nd real*8 w,u((2*nd+1)*ia,4,3),v(4),ar(3,3),ai(3,3),ww,t,gau complex*16 a(3,3) do 5 ix=1,3 do 5 jx=1,3 ar(ix,jx)=0.0d0 ai(ix,jx)=0.0d0 do 5 iv=is+1,is+2 ww=v(iv)**2-w**2 t=2.0d0*v(iv)*u(ik,iv,ix)*u(ik,iv,jx)/(ww**2+(v(iv)*gau)**2) ar(ix,jx)=ar(ix,jx) + ww * t 5 ai(ix,jx)=ai(ix,jx) + w*gau * t c wv^2-w^2 + iwG c a_ij = 2 wy ----------------------- ui uj c (wv^2-w^2)^2 + wv^2G^2) do 1 ix=1,3 do 1 jx=1,3 1 a(ix,jx)=dcmplx(ar(ix,jx),ai(ix,jx)) return end subroutine nv(v) real*8 v(*),a,sp a=dsqrt(sp(v,v)) v(1)=v(1)/a v(2)=v(2)/a v(3)=v(3)/a return end subroutine rda(ia,coa) implicit none integer*4 ia,coa(ia,6),i,j open(9,file='AM.SCR') do 1 i=1,ia 1 read(9,608)(coa(i,j),j=1,6) 608 format(5i6) close(9) return end subroutine raa(ia,r,q,nat) c r in Angstroms implicit none integer*4 ia,q(*),nat,i,j,io,ic,nn,k,jc(6),icn,ihn real*8 r(*),xi,yi,zi,xj,yj,zj,dij,xk,yk,zk,dkj ia=0 open(9,file='AM.SCR') do 1 i=1,nat xi=r(1+3*(i-1)) yi=r(2+3*(i-1)) zi=r(3+3*(i-1)) do 5 j=1,6 5 jc(j)=0 if(q(i).eq.6)then jc(1)=i io=0 ic=0 nn=0 icn=0 ihn=0 do 2 j=1,nat if(j.ne.i)then xj=r(1+3*(j-1)) yj=r(2+3*(j-1)) zj=r(3+3*(j-1)) dij=dsqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) if(dij.lt.1.65d0.and.q(j).eq.6)then ic=ic+1 jc(2)=j endif if(dij.lt.1.35d0.and.q(j).eq.8)then io=io+1 jc(3)=j endif if(dij.lt.1.55d0.and.q(j).eq.7)then nn=nn+1 jc(4)=j do 3 k=1,nat if(k.ne.i)then xk=r(1+3*(k-1)) yk=r(2+3*(k-1)) zk=r(3+3*(k-1)) dkj=dsqrt((xk-xj)**2+(yk-yj)**2+(zk-zj)**2) if(dkj.lt.1.11d0.and.q(k).eq.1)then ihn=ihn+1 jc(5)=k endif if(dkj.lt.1.55d0.and.q(k).eq.6)then icn=icn+1 jc(6)=k endif endif 3 continue endif endif 2 continue if(io.eq.1.and.ic.eq.1.and.nn.eq.1)then ia=ia+1 write(6,607)ia,(jc(j),j=1,6) 607 format(i6,':',i6,' CO',i6,' C',i6,' O',i6,' N',i6,' H',i6,' CN') c 3O c || c 1C 6C c / \ / c 2C 4N c | c 5H write(9,608)ia,(jc(j),j=1,6) 608 format(5i6) endif endif 1 continue close(9) return end subroutine getdips(N,nat,natu,nu3,uu,m1,r,iafirst,rq,ecm, 1lsmat,wmin,wmax,qvmi,qvma,itt,it,umin,umax,am,v, 1DSX,QSX,MSX,ds,rs,P,qq,A,Sqr,Sqi,qmax,delq,iq,ldip,AlPHA,G,AA, 1lttt,ram,ENM,iwr) implicit none integer*4 nat,natu,nu3,uu(nat,nat),ia,at,m1,ix,jx,iu,ii, 1iafirst(*),I,N,jj,it,iq,iy,iz,itt,umin,umax,j,kx,lx,iwr real*8 r(*),BOHR,ur(3),ui(3),qr(6),qi(6),mr(3),mi(3),ENM, 1sI,cI,pi,rq,re,im,FD,FC,debye,CM,AMU,clight,WR,ecm,auk, 1aulambda,tpn,qqr(3,3),qqi(3,3),wmin,wmax,qvmi,qvma, 1am(3,3),v(3),DSX(*),MSX(*),QSX(*),ds,rs,sp,ram(*),SFAC, 1P(nat,3,3),qq(3*nat,6),A(nat,3,3),Sqr(nu3,nu3),Sqi(nu3,nu3), 1qmax,delq,F,ALPHA(3*nat,3,3),G(3*nat,3,3),AA(3*nat,3,3,3), 1alr(3,3),ali(3,3),gr(3,3),gi(3,3),aar(3,3,3),aai(3,3,3), 1u1r(3),m1r(3),u1i(3),m1i(3) real*8,allocatable::Xcell(:),Pat(:,:),Aat(:,:), 1Qat(:,:),svr(:),svi(:),Acom(:,:),Qcom(:,:), 1AAcom(:,:,:,:),Gcom(:,:,:),ALat(:,:,:),Gat(:,:,:), 1AAat(:,:,:,:) logical lsmat,ldip,lttt allocate(Xcell(nu3),Pat(nu3,3),Aat(nu3,3),Qat(nu3,6), 1svr(nu3),svi(nu3),Acom(nu3,3),Qcom(nu3,6)) if(lttt)allocate(AAcom(nu3,3,3,3),Gcom(nu3,3,3)) SFAC=0.0234280d0 debye=2.541732d0 CM=219474.63d0 AMU=1822.887D0 clight=137.03604d0 FD =debye*dsqrt(CM/AMU/2.0d0) FC=debye*dsqrt(2.0d0/CM/AMU)/clight BOHR=0.52917705993d0 pi=4.0d0*atan(1.0d0) WR=dsqrt(ecm) tpn=2.0d0*pi*rq/dble(N) c wavelength in atomic units aulambda=1.0d0/ecm*1.0d8/BOHR c wave vector in atomic units auk=2.0d0*pi/aulambda c save atomic AATs to Aat,Atomic positions to xcell,P to Pat c in DOG coordinates: do 15 ia=1,natu at=uu(m1,ia) do 15 ix=1,3 ii=3*(ia-1)+ix Xcell(ii)=r(ix+3*(at-1))/BOHR do 15 jx=1,3 Pat(ii,jx )=P( at,ix,jx) Qat(ii,jx )=qq(ix+3*(at-1),jx ) Qat(ii,jx+3)=qq(ix+3*(at-1),jx+3) 15 Aat(ii,jx )=A( at,ix,jx) if(lttt)then c local atomic polarizability derivatives: allocate(ALat(nu3,3,3),Gat(nu3,3,3),AAat(nu3,3,3,3)) do 1 ia=1,natu at=uu(m1,ia) do 1 ix=1,3 ii=3*(ia-1)+ix do 1 jx=1,3 do 1 kx=1,3 ALat(ii,jx,kx)=ALPHA(ix+3*(at-1),jx,kx) Gat( ii,jx,kx)= G(ix+3*(at-1),jx,kx) do 1 lx=1,3 1 AAat(ii,jx,kx,lx)=AA(ix+3*(at-1),jx,kx,lx) endif c S-matrix vector for mode iq: do 19 ia=1,natu do 19 ix=1,3 iu=ix+3*(uu(m1,ia)-iafirst(m1)) ii=ix+3*(ia-1) svr(ii)=Sqr(iu,iq) 19 svi(ii)=Sqi(iu,iq) c zero dipole/polarizability moments: call vzd(ur,ui,qr,qi,mr,mi) if(lttt)call tz(alr,ali,gr,gi,aar,aai) c Loop over cell contributions: do 14 I=0,N-1 c phase factor: sI=sin(tpn*dble(I)) cI=cos(tpn*dble(I)) c put AAT in common system of I: call putaat(Xcell,Pat,Aat,Acom,natu) c put polarizability derivatives to sommon system: if(lttt)then if(iwr.eq.2)write(6,*)'Cell ',I,' local:' if(iwr.eq.2)call writetttd(natu,ALat,Gat,AAat) call TTT2(Xcell,ALat,Gat,AAat,Gcom,AAcom,natu) if(iwr.eq.2)call wrx(Xcell,natu) if(iwr.eq.2)write(6,*)'Cell ',I,' common:' if(iwr.eq.2)call writetttd(natu,ALat,Gcom,AAcom) endif c put Quadr. dervs. in common system of I: call putqat(Xcell,Pat,Qat,Qcom,natu) c contribution to the dipole moments: u1r=0.0d0 u1i=0.0d0 m1r=0.0d0 m1i=0.0d0 do 143 jj=1,nu3 re=svr(jj)*cI-svi(jj)*sI im=svr(jj)*sI+svi(jj)*cI if(lttt)then do 144 ix=1,3 do 144 jx=1,3 alr(ix,jx)=alr(ix,jx)+ALat(jj,ix,jx)*re ali(ix,jx)=ali(ix,jx)+ALat(jj,ix,jx)*im gr( ix,jx)=gr( ix,jx)+Gcom(jj,ix,jx)*re gi( ix,jx)=gi( ix,jx)+Gcom(jj,ix,jx)*im do 144 kx=1,3 aar(ix,jx,kx)=aar(ix,jx,kx)+AAcom(jj,ix,jx,kx)*re 144 aai(ix,jx,kx)=aai(ix,jx,kx)+AAcom(jj,ix,jx,kx)*im endif do 143 ix=1,3 qr(ix )=qr(ix )+Qcom(jj,ix )*re/WR*FD*0.5d0 qr(ix+3)=qr(ix+3)+Qcom(jj,ix+3)*re/WR*FD*0.5d0 qi(ix )=qi(ix )+Qcom(jj,ix )*im/WR*FD*0.5d0 qi(ix+3)=qi(ix+3)+Qcom(jj,ix+3)*im/WR*FD*0.5d0 m1r(ix)=m1r(ix)+Acom(jj,ix)*re m1i(ix)=m1i(ix)+Acom(jj,ix)*im u1r(ix)=u1r(ix)+Pat( jj,ix)*re u1i(ix)=u1i(ix)+Pat( jj,ix)*im mr(ix)=mr(ix)+Acom(jj,ix)*re*WR*FC mi(ix)=mi(ix)+Acom(jj,ix)*im*WR*FC ur(ix)=ur(ix)+Pat( jj,ix)*re/WR*FD 143 ui(ix)=ui(ix)+Pat( jj,ix)*im/WR*FD call rewritew(qqr,qr,qqi,qi) c S-matrix,temporary record: call wrs(natu,lsmat,ecm,rq,wmin,wmax,qvmi,qvma,itt,it,I,svr,svi, 1umin,umax,sI,cI) c prepare for the next: c rotate local tensors forward call trr(Pat,natu,am) call trr(Aat,natu,am) call trq(Qat,natu,am) call trs(natu,svr,am) call trs(natu,svi,am) if(lttt)then call tra(ALat,natu,am) call tra(Gat ,natu,am) call trt(AAat,natu,am) endif c rotate and advance coordinates (Xcell in bohrs and v in A): 14 call trc(Xcell,am,natu,v) if(lttt)then c multiply by SFAC = 1 /sqrt(AMU) for compatibility with new6: do 147 ix=1,3 do 147 jx=1,3 alr(ix,jx)=alr(ix,jx)*SFAC ali(ix,jx)=ali(ix,jx)*SFAC gr( ix,jx)=gr( ix,jx)*SFAC gi( ix,jx)=gi( ix,jx)*SFAC do 147 kx=1,3 aar(ix,jx,kx)=aar(ix,jx,kx)*SFAC 147 aai(ix,jx,kx)=aai(ix,jx,kx)*SFAC if(iwr.ge.2)call wrl(iq,ecm,alr,ali,gr,gi,aar,aai) endif if(ldip)write(29,619)iq,ecm,ur,ui,(1000.0d0*mr(ix),ix=1,3), 1(1000.0d0*mi(ix),ix=1,3), 1((qqr(i,j),i=1,j),j=1,3), 1((qqi(i,j),i=1,j),j=1,3) 619 format(i3,f8.1,' u:',6f10.4,/,11x,' m:',6f10.4,/, 111x,' q:',6f10.4,/,11x,' q:',6f10.4) c f=delq/qmax c for X, Y and Z beam directions: do 58 ix=1,3 iy=ix+1 if(iy.eq.4)iy=1 iz=iy+1 if(iz.eq.4)iz=1 DSX(ix)= 1.5d0*(ur(iy)*ur(iy) +ur(iz)*ur(iz) 1 +ui(iy)*ui(iy) +ui(iz)*ui(iz) )*F QSX(ix)=-auk*1.5d0*(ur(iy)*qqr(iz,ix)-ur(iz)*qqr(ix,iy) 1 +ui(iy)*qqi(iz,ix)-ui(iz)*qqi(ix,iy))*F 58 MSX(ix)= 1.5d0*(ur(iy)*mr(iy) +ur(iz)*mr(iz) 1 +ui(iy)*mi(iy) +ui(iz)*mi(iz) )*F c dipole strength: ds=(sp(ur,ur)+sp(ui,ui))*F c rotational strength: rs=(sp(ur,mr)+sp(ui,mi))*F c c raman: if(lttt)call groa(ram,alr,ali,gr,gi,aar,aai,ENM,iwr) return end subroutine getdipr(iqt,nu3,prc,pic,arc,aic,qrc,qic,Sqr,Sqi,N,ecm, 1rq,tau,ds,rs,DSX,MSX,QSX,v,umr,umi,at,ldip,lttt,alrc,alic,grc, 1gic,aarc,aaic,ENM,ram,iwr) implicit none integer*4 nu3,iq,N,i,j,k,iqt,l,iwr real*8 prc(nu3,3),pic(nu3,3),arc(nu3,3),aic(nu3,3),rq,sr,si, 1qrc(nu3,3,3),qic(nu3,3,3),Sqr(nu3,nu3),Sqi(nu3,nu3),tau,l2,ur(3), 1ui(3),mr(3),mi(3),qr(3,3),qi(3,3),pi,AN,urt(3),uit(3), 1mrt(3),mit(3),qrt(3,3),qit(3,3),ak,al,ecm,debye,CM,AMU,clight,FD, 1FC,BOHR,WR,aulambda,auk,ds,rs,DSX(*),MSX(*),QSX(*),s2r,s2i, 1v(3),sp,f1,urc(3),uic(3),qro(3,3),qio(3,3),ENM,ram(*),l1, 1mrc(3),mic(3),umr(3,3),umi(3,3),at(3,3),alrc(nu3,3,3),aj, 1alic(nu3,3,3),grc(nu3,3,3),gic(nu3,3,3),aarc(nu3,3,3,3), 1aaic(nu3,3,3,3),ali(3,3),alr(3,3),gi(3,3),gr(3,3),aai(3,3,3), 1aar(3,3,3),alti(3,3),altr(3,3),gti(3,3),gtr(3,3),aati(3,3,3), 1aatr(3,3,3),alio(3,3),alro(3,3),gio(3,3),gro(3,3),aaio(3,3,3), 1aaro(3,3,3),ll logical ldip,lttt iq=iqt AN=dble(N) pi=4.0d0*datan(1.0d0) f1=2.0d0*pi*rq/AN debye=2.541732d0 CM=219474.63d0 AMU=1822.887D0 clight=137.03604d0 FD =debye*dsqrt(CM/AMU/2.0d0) FC=debye*dsqrt(2.0d0/CM/AMU)/clight BOHR=0.52917705993d0 al=dsqrt(sp(v,v))/4.0d0/BOHR l2=al*2.0d0 l1=l2*2.0d0 WR=dsqrt(ecm) c wavelength in atomic units aulambda=1.0d0/ecm*1.0d8/BOHR c wave vector in atomic units auk=2.0d0*pi/aulambda c transform complex dipolar derivatives to normal mode coordinates c du/dQ = dU/dx . S call trpa(iq,nu3,prc,pic,Sqr,Sqi,ur,ui) call trpa(iq,nu3,arc,aic,Sqr,Sqi,mr,mi) call trpq(iq,nu3,qrc,qic,Sqr,Sqi,qr,qi) if(lttt)then call trpq(iq,nu3,alrc,alic,Sqr,Sqi,alr,ali) call trpq(iq,nu3,grc,gic,Sqr,Sqi,gr,gi) call traa(iq,nu3,aarc,aaic,Sqr,Sqi,aar,aai) endif do 1 k=1,3 ak=dble(k-2) call s1s(f1+ak*tau,sr,si,AN) urt(k)=(ur(k)*sr-ui(k)*si)/WR*FD uit(k)=(ui(k)*sr+ur(k)*si)/WR*FD mrt(k)=(mr(k)*sr-mi(k)*si)*WR*FC mit(k)=(mi(k)*sr+mr(k)*si)*WR*FC do 1 j=1,3 aj=dble(j-2) call s1s(f1+(ak+aj)*tau,sr,si,AN) qrt(k,j)=(qr(k,j)*sr-qi(k,j)*si)/WR*FD*0.5d0 1 qit(k,j)=(qi(k,j)*sr+qr(k,j)*si)/WR*FD*0.5d0 if(lttt)then do 11 k=1,3 ak=dble(k-2) do 11 j=1,3 aj=dble(j-2) call s1s(f1+(ak+aj)*tau,sr,si,AN) altr(k,j)=(alr(k,j)*sr-ali(k,j)*si) alti(k,j)=(ali(k,j)*sr+alr(k,j)*si) gtr( k,j)=(gr( k,j)*sr-gi( k,j)*si) gti( k,j)=(gi( k,j)*sr+gr( k,j)*si) do 11 l=1,3 ll=dble(l-2) call s1s(f1+(ak+aj+ll)*tau,sr,si,AN) aatr(k,j,l)=(aar(k,j,l)*sr-aai(k,j,l)*si) 11 aati(k,j,l)=(aai(k,j,l)*sr+aar(k,j,l)*si) endif c add non-local parts c m( 1): 1 0 1 mi = i eps_-ijk Rj uk call s2s(f1+tau,s2r,s2i,AN) mrt(3)=mrt(3)-al*(ui(3)*s2r+ur(3)*s2i)*WR*FC mit(3)=mit(3)+al*(ur(3)*s2r-ui(3)*s2i)*WR*FC c m(-1): -1 0 -1 call s2s(f1-tau,s2r,s2i,AN) mrt(1)=mrt(1)+al*(ui(1)*s2r+ur(1)*s2i)*WR*FC mit(1)=mit(1)-al*(ur(1)*s2r-ui(1)*s2i)*WR*FC c 00 = 2 u0 L S2(0): call s2s(f1,s2r,s2i,AN) qrt(2,2)=qrt(2,2)+2.0d0*l2*(ur(2)*s2r-ui(2)*s2i)/WR*FD qit(2,2)=qit(2,2)+2.0d0*l2*(ui(2)*s2r+ur(2)*s2i)/WR*FD c -1 0 = u_-1 L S2(-1) call s2s(f1-tau,s2r,s2i,AN) qrt(1,2)=qrt(1,2)+ l2*(ur(1)*s2r-ui(1)*s2i)/WR*FD qit(1,2)=qit(1,2)+ l2*(ui(1)*s2r+ur(1)*s2i)/WR*FD c 1 0 = u_ 1 L S2( 1) call s2s(f1+tau,s2r,s2i,AN) qrt(3,2)=qrt(3,2)+ l2*(ur(3)*s2r-ui(3)*s2i)/WR*FD qit(3,2)=qit(3,2)+ l2*(ui(3)*s2r+ur(3)*s2i)/WR*FD qrt(2,1)=qrt(1,2) qit(2,1)=qit(1,2) qrt(2,3)=qrt(3,2) qit(2,3)=qit(3,2) if(lttt)then do 2 i=1,3 c i -1 = i x (s(i-1) alpha_i -1) call s2s(f1+dble(i-2-1)*tau,s2r,s2i,AN) gtr(i,1)=gtr(i,1) - l2*(alr(i,1)*s2i+ali(i,1)*s2r) gti(i,1)=gti(i,1) + l2*(alr(i,1)*s2r-ali(i,1)*s2i) c i 1 =-i x (s(i 1) alpha_i 1) call s2s(f1+dble(i-2+1)*tau,s2r,s2i,AN) gtr(i,3)=gtr(i,3) + l2*(alr(i,3)*s2i+ali(i,3)*s2r) 2 gti(i,3)=gti(i,3) - l2*(alr(i,3)*s2r-ali(i,3)*s2i) do 3 i=1,3 do 3 j=1,3 do 3 k=1,3 c delta(aj,-ak): if(j+k.eq.4)then call s2s(f1+dble(i-2)*tau,s2r,s2i,AN) aatr(i,j,k)=aatr(i,j,k) - l1*(alr(i,2)*s2r-ali(i,2)*s2i) aati(i,j,k)=aati(i,j,k) - l1*(ali(i,2)*s2r+alr(i,2)*s2i) endif if(j.eq.2)then call s2s(f1+dble(i-2+k-2)*tau,s2r,s2i,AN) aatr(i,j,k)=aatr(i,j,k) + 1.5d0*l1*(alr(i,k)*s2r-ali(i,k)*s2i) aati(i,j,k)=aati(i,j,k) + 1.5d0*l1*(ali(i,k)*s2r+alr(i,k)*s2i) endif if(k.eq.2)then call s2s(f1+dble(i-2+j-2)*tau,s2r,s2i,AN) aatr(i,j,k)=aatr(i,j,k) + 1.5d0*l1*(alr(i,j)*s2r-ali(i,j)*s2i) aati(i,j,k)=aati(i,j,k) + 1.5d0*l1*(ali(i,j)*s2r+alr(i,j)*s2i) endif 3 continue endif c transform complex dipole to new cartesians ut -> uc call toc(urc,uic,urt,uit,umr,umi) call toc(mrc,mic,mrt,mit,umr,umi) call toQ(qro,qio,qrt,qit,umr,umi) if(lttt)then call toQ(alro,alio,altr,alti,umr,umi) call toQ( gro, gio, gtr, gti,umr,umi) call to3(aaro,aaio,aatr,aati,umr,umi) endif c transform dipole new cartesians to old cartesians uc -> uc call tto(urc,uic,at) call tto(mrc,mic,at) call ttQ(qro,qio,at) if(lttt)then call ttQ(alro,alio,at) call ttQ( gro, gio,at) call tt3(aaro,aaio,at) endif if(ldip)write(29,619)iq,ecm, 1(urc(i),i=1,3),(uic(i),i=1,3), 1(mrc(i),i=1,3),(mic(i),i=1,3), 1((qro(i,j)/10.0d0,i=1,j),j=1,3), 1((qio(i,j)/10.0d0,i=1,j),j=1,3) 619 format(i3,f8.1,' u:',6f10.4,/,11x,' m:',6f10.4,/, 111x,' q:',6f10.4,/,11x,' q:',6f10.4) c q: xx xy yy xz yz zz ds=sp(urc,urc)+sp(uic,uic) rs=sp(urc,mrc)+sp(uic,mic) c for i, j and k beam directions: do 58 i=1,3 j=i+1 if(j.eq.4)j=1 k=j+1 if(k.eq.4)k=1 DSX(i)= 1.5d0*(urc(j)*urc(j) + urc(k)*urc(k) 1 +uic(j)*uic(j) + uic(k)*uic(k) ) QSX(i)=-auk*1.5d0*(urc(j)*qro(k,i)-urc(k)*qro(i,j) 1 +uic(j)*qio(k,i)-uic(k)*qio(i,j)) 58 MSX(i)= 1.5d0*(urc(j)*mrc(j) + urc(k)*mrc(k) 1 +uic(j)*mic(j) + uic(k)*mic(k) ) if(lttt)call groa(ram,alro,alio,gro,gio,aaro,aaio,ENM,iwr) return end subroutine groa(ram,alr,ali,gr,gi,ar,ai,ENM,iwr) implicit none integer*4 I,J,K,L,iwr real*8 alr(3,3),ali(3,3),gr(3,3),gi(3,3),ar(3,3,3),ai(3,3,3), 1ram(*),SAL0,SAL1,SAG0,SAG1,SA1,SPALr,SPALi,EPS, 1S180,D180,D,YDX,YDY,P1,doc,AMU,BOHR,SpA3r,SpA3i,SpA32, 1gpisvejc,beta2,betaa,betag,roa3,roa1,ram1, 1EXCA,ENM EXCA=ENM*10.0d0 AMU=1822.0d0 BOHR=0.529177d0 c gpisvejc=980.0 gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 SPALr=0.0d0 SPALi=0.0d0 DO 3 I=1,3 SPALr=SPALr+alr(I,I) SPALi=SPALi+ali(I,I) DO 3 J=1,3 SAL0=SAL0+alr(I,I)*alr(J,J)+ali(I,I)*ali(J,J) SAL1=SAL1+alr(I,J)*alr(I,J)+ali(I,J)*ali(I,J) SAG0=SAG0+alr(I,I)*gr( J,J)+ali(I,I)*gi( J,J) SAG1=SAG1+alr(I,J)*gr( I,J)+ali(I,J)*gi( I,J) DO 3 K=1,3 DO 3 L=1,3 3 SA1=SA1+EPS(I,K,L)*(alr(I,J)*ar(K,L,J)+ali(I,J)*ai(K,L,J))/3.0d0 if(iwr.ge.2)write(6,609)SAL0,SAL1,SAG0,SAG1,SA1 609 format(' SAL0 SAL1 SAG0 SAG1 SA1:',5g11.4) c IL+/-IR, backscattering S180=7.0d0*SAL1+SAL0 c write(6,*)ECM,SAL0,SAL1,SAG0,SAG1,SA1 D180=8.0d0*(3.0d0*SAG1-SAG0+SA1) D=D180 IF(S180.gt.1.0d-9)D=D/S180 c c polarized backscattering components (YDY+YDX=S180): YDY=6.0d0*(3.0d0*SAL1- SAL0)*(AMU*BOHR**4) YDX=6.0d0*(4.0d0*SAL1+2.0d0*SAL0)*(AMU*BOHR**4) if(YDX.gt.1.0d-9)then P1=YDY/YDX else P1=0.0d0 endif SpA3r=SPALr/3.0d0*dsqrt((AMU*BOHR**4)) SpA3i=SPALi/3.0d0*dsqrt((AMU*BOHR**4)) SpA32=SpA3r**2+SpA3i**2 beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) c c Barron's tensor invariants: c alpha G = [(1/3)Sp alpha] [(1/3)Sp G'] etc c beta(A)^2=Delta2 in Gaussian: betaa=SA1*3.0d0/2.0d0 *gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0 *gpisvejc c c DCP 180 c ram3=24.0d0*beta2 roa3=8.0d0*(12.0d0*betag+4.0d0*betaa) c c ICP(180): c 8/c*(12*beta(G')^2 + 4*beta(A)^2): roa1=96.0d0*(betag+betaa/3.0d0) ram1=4.0d0*(45.0d0*SpA32+7.0d0*beta2) doc=(3*P1-1)/(P1+1) c Raman, parallel polarization: ram(1)=YDX c Raman, perpendicular poalrization: ram(2)=YDY c Raman, total: ram(3)=ram1 c Depolarization ratio: ram(4)=P1 c ROA SCP 180: ram(5)=roa1 c DOC: ram(6)=doc c DCP ROA 180: ram(7)=roa3 do 6 i=1,3 call setjk(i,j,k) c S0(x)=(1/2)(azz azz + ayy ayy + azy azy + azy azy) c S0(y)=(1/2)(axx axx + azz azz + axz axz + axz axz) c S0(z)=(1/2)(axx axx + ayy ayy + axy axy + axy axy) c jj jj kk kk jk jk jk jk c 8-10 oriented ram1 ram(7+i)=180.0d0*0.5d0*AMU*BOHR**4*( 1alr(j,j)*alr(j,j)+alr(k,k)*alr(k,k)+ 1alr(j,k)*alr(j,k)+alr(j,k)*alr(j,k)+ 1ali(j,j)*ali(j,j)+ali(k,k)*ali(k,k)+ 1ali(j,k)*ali(j,k)+ali(j,k)*ali(j,k)) c -2ayx(Gxy+Gyx)+(ayy-axx)(Gxx-Gyy) c kj jk kj kk jj jj kk c +(2w/3)(ayyAx zy-axxAy zx+ayxAx zx-axyAy zy) c kk j ik jj k ij kj j ij jk k ik c 11-13 oriented roa1 6 ram(10+i)=-180.0d0*gpisvejc*( 1-2.0d0*alr(j,k)*(gr(j,k)+gr(k,j)) 1-2.0d0*ali(j,k)*(gi(j,k)+gi(k,j)) 1+(alr(k,k)-alr(j,j))*(gr(j,j)-gr(k,k)) 1+(ali(k,k)-ali(j,j))*(gi(j,j)-gi(k,k)) 1+2.0d0/3.0d0* 1(alr(k,k)*Ar(j,i,k)-alr(j,j)*Ar(k,i,j) 1+alr(k,j)*Ar(j,i,j)-alr(j,k)*Ar(k,i,k) 1+ali(k,k)*Ai(j,i,k)-ali(j,j)*Ai(k,i,j) 1+ali(k,j)*Ai(j,i,j)-ali(j,k)*Ai(k,i,k))) RETURN END subroutine spreadr(ecm,ram,s,np,wsmin,wsmax,fwhm,lg,temp) implicit none real*8 s(16,np),wsmin,wsmax,w,dw,ecm,pd,fwhm,sf,ram(*), 1dd(8),prod,temp,c integer*4 i,np,j logical lg c=1.0d0/ecm/(1.0d0-exp(-(1.44d0*ecm/temp))) dd(1)=ram( 3)*c dd(2)=ram( 8)*c dd(3)=ram( 9)*c dd(4)=ram(10)*c dd(5)=ram( 5)*c dd(6)=ram(11)*c dd(7)=ram(12)*c dd(8)=ram(13)*c dw=(wsmax-wsmin)/dble(np-1) w=wsmin-dw if(lg)then pd=fwhm/2.0d0/dsqrt(log(2.0d0)) else pd=fwhm/2.0d0 endif do 1 i=1,np w=w+dw prod=sf(pd,lg,w,ecm,1) do 1 j=1,8 1 s(8+j,i)=s(8+j,i)+prod*dd(j) return end subroutine sspread(ecm,d1,d2,d3,d4,d5,d6,d7,d8, 1s,np,wsmin,wsmax,fwhm,lg) implicit none real*8 s(16,np),wsmin,wsmax,w,dw,ecm,pd,fwhm,sf,d1, 1d2,d3,d4,d5,d6,d7,d8,dd(8),prod integer*4 i,np,j logical lg dd(1)=d1*108.7d0 dd(2)=d2*108.7d0 dd(3)=d3*108.7d0 dd(4)=d4*108.7d0 dd(5)=d5*435.0d0 dd(6)=d6*435.0d0 dd(7)=d7*435.0d0 dd(8)=d8*435.0d0 dw=(wsmax-wsmin)/dble(np-1) w=wsmin-dw if(lg)then pd=fwhm/2.0d0/dsqrt(log(2.0d0)) else pd=fwhm/2.0d0 endif do 1 i=1,np w=w+dw prod=w*sf(pd,lg,w,ecm,1) do 1 j=1,8 1 s(j,i)=s(j,i)+prod*dd(j) return end subroutine wspectra(f,wsmin,wsmax,s,np,j) implicit none real*8 s(16,np),wsmin,wsmax,w,dw integer*4 i,np,j character*(*) f OPEN(30,FILE=f) dw=(wsmax-wsmin)/dble(np-1) w=wsmin-dw do 1 i=1,np w=w+dw 1 write(30,300)w,s(j,i) 300 format(f12.2,' ',g12.4) close(30) write(6,*)f//' written' return end subroutine smatout(ntt,umin,umax,natu,r,N,uu,nat,m1,am,v, 1nu3,z) implicit none integer*4 ntt,umin,umax,nut,I,ia,ix,natu,N,uu(nat,nat), 1nat,m1,NQ,nu3,im,jdum,i1,iu,z(*) real*8 r(*),am(3,3),v(*),ecm,BOHR real*8,allocatable::Xcell(:),smat(:,:),w(:) BOHR=0.52917705993d0 nut=umax-umin+1 NQ=ntt/nut allocate(w(NQ),Xcell(3*natu),smat(3*natu*nut,NQ)) write(6,607)NQ 607 format(' Writing S-matrix,',i7,' modes') rewind 46 open(44,file='FR.INP') write(44,400)NQ,3*natu*nut,natu*nut 400 FORMAT(3I7) do 151 ia=1,natu do 151 ix=1,3 151 Xcell(ix+3*(ia-1))=r(ix+3*(uu(m1,ia)-1))/BOHR do 142 I=0,N-1 if(I+1.ge.umin.and.I+1.le.umax)then do 152 ia=1,natu 152 write(44,440)z(uu(m1,ia)),(Xcell(ix+3*(ia-1))*BOHR,ix=1,3) 440 format(i7,3f12.6) endif c (Xcell in bohrs and v in A): 142 call trc(Xcell,am,natu,v) WRITE(44,140) 140 FORMAT(' Atom Mode X-disp. Y-disp. Z-disp.') do 141 ia=1,nu3*nut do 141 ix=1,NQ 141 smat(ia,ix)=0.0d0 c first read cell: i1=-99 im=0 1000 read(46,*,err=66,end=66)ecm,I,ecm,ecm if(i1.eq.-99)i1=I c record mode for the first read cell: if(I.eq.i1)then im=im+1 w(im)=ecm endif do 201 ia=1,natu do 201 ix=1,3 201 read(46,*)(smat(ix+3*(ia-1)+nu3*I,im),jdum=1,3) goto 1000 66 close(46) DO 22 I=umin-1,umax-1 DO 22 ia=1,natu DO 22 im=1,NQ IU=3*(ia-1)+nu3*I 22 WRITE(44,45)ia+natu*I,im,(smat(IU+ix,im),ix=1,3) 45 FORMAT(2I7,3F11.6) WRITE(44,*)NQ do 4 im=1,NQ WRITE(44,134)w(im) 134 format(F11.3,$) 4 if(mod(im,6).eq.0)write(44,*) close(44) return end subroutine mkfirst(nat,nu,iafirst,u,natu) IMPLICIT none integer*4 nat,nu,iafirst(*),u(nat,nat),i,j,natu c first atom in each init: do 3 i=1,nu iafirst(i)=u(i,1) do 3 j=2,natu 3 if(u(i,j).lt.iafirst(i))iafirst(i)=u(i,j) return end subroutine addtdc(q,nat,natu,n3,nu3,nu,f,m1,Dc,Ds,N,u, 1ru,A,v,P,iafirst,epsr,nd) IMPLICIT none real*8 f(n3,n3),Dc(nu3,nu3),Ds(nu3,nu3),pi,sJ,cJ,q,f1,rj(3), 1A(3,3),v(3),P(nat,3,3),epsr,bohr,xj,yj,zj,xi,yi,zi,d1,d2,d5, 1pix,piy,piz,pjx,pjy,pjz,ri(3),uij,uir,ujr,did,ru(*),pj(3,3) integer*4 nu,m1,nu3,j,nat,natu,u(nat,nat),jx,iu,ju,jb,ib, 1ja,N,N3,ia,ix,iafirst(*),nd pi=4.0d0*datan(1.0d0) f1=2.0d0*pi*q/dble(N) bohr=0.529177d0 write(6,*)' Adding TDC' c loop over -nd..nd elementary cells around m1 do 1 J=m1-nd,m1+nd c phase relative to m1: sJ=sin(f1*dble(J-m1)) cJ=cos(f1*dble(J-m1)) c exp( i 2 pi q I / N )/sqrt(N) on I = 0 ... N-1 do 1 ja=1,natu call asx(rj,ru,ja) c rotate to cell J: call rtcj(rj,A,J-m1) c position in cell J (relative to T1) xj=rj(1)+v(1)*dble(J-m1) yj=rj(2)+v(2)*dble(J-m1) zj=rj(3)+v(3)*dble(J-m1) do 1 ia=1,ja-1 call asx(ri,ru,ia) xi=(xj-ri(1))/bohr yi=(yj-ri(2))/bohr zi=(zj-ri(3))/bohr d2=xi**2+yi**2+zi**2 d1=dsqrt(d2) d5=d2*d2*d1 do 1 ix=1,3 c index in FILE.X ib=ix+3*(u(m1,ia)-1) c relative to the first atom in the cell iu=ix+3*(u(m1,ia)-iafirst(m1)) pix=p(u(m1,ia),ix,1) piy=p(u(m1,ia),ix,2) piz=p(u(m1,ia),ix,3) do 2 jx=1,3 pj(jx,1)=p(u(m1,ja),jx,1) pj(jx,2)=p(u(m1,ja),jx,2) 2 pj(jx,3)=p(u(m1,ja),jx,3) call rtc2(pj,A,J-m1) do 1 jx=1,3 pjx=pj(jx,1) pjy=pj(jx,2) pjz=pj(jx,3) ju=jx+3*(u(m1,ja)-iafirst(m1)) uij=pix*pjx+piy*pjy+piz*pjz uir=pix*xi +piy*yi +piz*zi ujr=pjx*xi +pjy*yi +pjz*zi did=(uij*d2-3.0d0*uir*ujr)/d5/epsr if(J.ge.1.and.J.le.nu)then c within supplied geometry jb=jx+3*(u(J,ja)-1) if(dabs(f(ib,jb)).lt.1.0d-7)then Dc(iu,ju)=Dc(iu,ju)+did*cJ Ds(iu,ju)=Ds(iu,ju)+did*sJ endif else Dc(iu,ju)=Dc(iu,ju)+did*cJ Ds(iu,ju)=Ds(iu,ju)+did*sJ endif 1 continue return end subroutine maked(q,nat,natu,n3,nu3,nu,f,m1,Dc,Ds,N,u,iwr,iafirst) IMPLICIT none real*8 f(n3,n3),Dc(nu3,nu3),Ds(nu3,nu3),pi,sJ,cJ,q,f1,mJ integer*4 nu,m1,nu3,j,nat,natu,u(nat,nat),jx,ii,jj,iu,ju, 1ja,N,N3,iwr,ia,ix,iafirst(*) pi=4.0d0*datan(1.0d0) f1=2.0d0*pi*q/dble(N) Dc=0.0d0 Ds=0.0d0 c loop over all available elementary cells: do 1 J=1,nu c phase relative to m1: mJ=f1*dble(J-m1) sJ=sin(mJ) cJ=cos(mJ) c exp( i 2 pi q I / N )/sqrt(N) on I = 0 ... N-1 do 1 ja=1,natu do 1 jx=1,3 c global index: jj=jx+3*(u( J,ja) -1) c relative to the first atom in the cell ju=jx+3*(u( J,ja)-iafirst(J)) do 1 ia=1,natu do 1 ix=1,3 ii=ix+3*(u(m1,ia) -1) iu=ix+3*(u(m1,ia)-iafirst(m1)) Dc(iu,ju)=Dc(iu,ju)+f(ii,jj)*cJ 1 Ds(iu,ju)=Ds(iu,ju)+f(ii,jj)*sJ if(iwr.gt.0)then write(6,602)Dc(1,1),Dc(2,2),Ds(1,1),Ds(2,2) 602 format(2f12.6,/,2f12.6) write(6,*)' D done' endif return end SUBROUTINE MASSFF(N,f,m) IMPLICIT none integer*4 i,j,N REAL*8 m(*),F(N,N) do 1 i=1,N do 1 j=1,N 1 f(i,j)=f(i,j)/dsqrt(m((i+2)/3)*m((j+2)/3)) WRITE(*,*)' Mass-scaled FF done.' return end SUBROUTINE READFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) OPEN(20,FILE='FILE.FC',STATUS='OLD') N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) CONST=4.359828d0/0.5291772d0/0.5291772d0 DO 3 I=1,N DO 3 J=I,N 3 FCAR(J,I)=FCAR(J,I)*CONST DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in.' CLOSE(20) RETURN END subroutine READM(nat,m) implicit none integer*4 i,nat real*8 m(*) open(9,file='FTRY.INP', status='old') read(9,*) read(9,*) read(9,90)(m(i),i=1,nat) 90 format(6F12.6) read(9,*) read(9,90)(m(i),i=1,nat) close(9) WRITE(*,*)' Masses read in.' return end subroutine crms(r,nat,natu,u,m1,m2,v,a,T1) IMPLICIT none integer*4 i,nat,natu,u(nat,nat),ix,J,m1,m2 real*8 T1(3),T2(3),v(3),r(*),ANG(3),a(3,3) real*8,allocatable::X1(:,:),X2(:,:) allocate(X1(natu,3),X2(natu,3)) c rewrite coordinates: do 1 i=1,natu do 1 ix=1,3 X1(i,ix)=r(ix+3*(u(m1,i)-1)) 1 X2(i,ix)=r(ix+3*(u(m2,i)-1)) C Calculate center of the coordinates that will be rotated: DO 33 J=1,3 T2(J) =0.0d0 T1(J)=0.0d0 DO 33 I=1,natu T2(J)=T2(J)+X2(I,J) 33 T1(J)=T1(J)+X1(I,J) DO 38 J=1,3 T2(J)=T2(J)/DBLE(natu) 38 T1(J)=T1(J)/DBLE(natu) C Subtract mass-center: DO 35 I=1,natu DO 35 J=1,3 X1(I,J)=X1(I,J)-T1(J) 35 X2(I,J)=X2(I,J)-T2(J) do 2 ix=1,3 2 V(ix)=T2(ix)-T1(ix) write(6,601)(v(ix),ix=1,3) 601 format(' Shift vector: ',/,3f10.4) A(1,1)=1.0d0 CALL DOMATRIX(natu,A,X1,X2,ANG) call wm('Transformation matrix:',A) write(6,603)(180.0d0*ANG(ix)/4.0d0/datan(1.0d0),ix=1,3) 603 format(' Angles:',/,3f10.4) c R2a = T1a + va + Aab (R1 - T1)b, etc. return end subroutine crmi(r,nat,natu,u,m1,m2,v,a,T1,tau) IMPLICIT none integer*4 i,nat,natu,u(nat,nat),ix,jx,J,m1,m2,ierr real*8 T1(3),T2(3),v(3),r(*),TH1(3,3),a(3,3), 1TH2(3,3),A1(3,3),A2(3,3),TI1(3),TI2(3),v1(3),v2(3), 1sp,pi,zv(3),vy(3),vx(3),tau,v12(3) real*8,allocatable::X1(:,:),X2(:,:) allocate(X1(natu,3),X2(natu,3)) pi=4.0d0*datan(1.0d0) c rewrite coordinates: do 1 i=1,natu do 1 ix=1,3 X1(i,ix)=r(ix+3*(u(m1,i)-1)) 1 X2(i,ix)=r(ix+3*(u(m2,i)-1)) C Calculate center of the coordinates that will be rotated: DO 33 J=1,3 T2(J) =0.0d0 T1(J)=0.0d0 DO 33 I=1,natu T2(J)=T2(J)+X2(I,J) 33 T1(J)=T1(J)+X1(I,J) DO 38 J=1,3 T2(J)=T2(J)/DBLE(natu) 38 T1(J)=T1(J)/DBLE(natu) do 2 ix=1,3 2 V(ix)=T2(ix)-T1(ix) write(6,601)(v(ix),ix=1,3) 601 format(' Shift vector: ',/,3f10.4) C Subtract mass-center: DO 35 I=1,natu DO 35 J=1,3 X1(I,J)=X1(I,J)-T1(J) 35 X2(I,J)=X2(I,J)-T2(J) c momenta of inertia do 3 ix=1,3 do 3 jx=1,3 TH1(ix,jx)=0.0d0 TH2(ix,jx)=0.0d0 do 3 i=1,natu TH1(ix,jx)=TH1(ix,jx)+X1(i,ix)*X1(i,jx) 3 TH2(ix,jx)=TH2(ix,jx)+X2(i,ix)*X2(i,jx) write(6,605)((TH1(ix,jx),ix=1,3),jx=1,3) 605 format(' Moment of inertia I :',/,3(3f10.4,/)) ierr=0 CALL TRED12(3,TH1,A1,TI1,2,IERR) if(ierr.ne.0)call report('cannot diagonalize') write(6,603)(TI1(ix),ix=1,3) 603 format(' Eigenvalues :',3f10.4) write(6,600)((TH2(ix,jx),ix=1,3),jx=1,3) 600 format(' Moment of inertia II:',/,3(3f10.4,/)) ierr=0 CALL TRED12(3,TH2,A2,TI2,2,IERR) if(ierr.ne.0)call report('cannot diagonalize') write(6,603)(TI2(ix),ix=1,3) c v1,v2 ... largest moment of inertia vector x shift v1(1)=A1(2,1)*V(3)-A1(3,1)*V(2) v1(2)=A1(3,1)*V(1)-A1(1,1)*V(3) v1(3)=A1(1,1)*V(2)-A1(2,1)*V(1) v2(1)=A2(2,1)*V(3)-A2(3,1)*V(2) v2(2)=A2(3,1)*V(1)-A2(1,1)*V(3) v2(3)=A2(1,1)*V(2)-A2(2,1)*V(1) call norm(v1) call norm(v2) tau=acos(sp(v1,v2)) call vp(v12,v1,v2) if(sp(v12,v).lt.0.0d0)tau=-tau write(6,604)tau*180.0d0/pi 604 format(' Twist angle:',f10.2,' deg') c try to use old axes as much as possible: c set new z-coordinate along v: zv=v call norm(zv) c nex x axis along old x or y, as far as they do not go along v: vx=0.0d0 if(dabs(v(1)).gt.dabs(v(2)).and.dabs(v(1)).gt.dabs(v(3)))then vx(2)=1.0d0 else vx(1)=1.0d0 endif call vp(vy,zv,vx) call norm(vy) call vp(vx,vy,zv) call norm(vx) c transformation matrix - first index old do 5 ix=1,3 A(ix,1)=vx(ix) A(ix,2)=vy(ix) 5 A(ix,3)=zv(ix) call wm('Trasformation matrix to new Cartesians:',a) return end subroutine vp(vy,vz,vx) real*8 vy(3),vz(3),vx(3) vy(1)=vz(2)*vx(3)-vz(3)*vx(2) vy(2)=vz(3)*vx(1)-vz(1)*vx(3) vy(3)=vz(1)*vx(2)-vz(2)*vx(1) return end SUBROUTINE DOMATRIX(natu,A,XS,XB,ANG) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) real*8 A(3,3),XS(natu,3),XB(natu,3) C STARTING VALUES FOR THE ITERATION: ANG(1)=2.1d0 ANG(2)=0.0d0 ANG(3)=0.0d0 RMS=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG,A,natu,XS,XB,RMS) ANG(1)=0.0d0 ANG(2)=0.1d0 ANG(3)=2.1d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG,A,natu,XS,XB,RMS) ANG(1)=0.2d0 ANG(2)=2.1d0 ANG(3)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG,A,natu,XS,XB,RMS) ANG(1)=1.4d0 ANG(2)=0.0d0 ANG(3)=0.1d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG,A,natu,XS,XB,RMS) FTOL=1.0d-9 ITMAX=1500 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-10)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)then RETURN endif IF(ITER.EQ.ITMAX)call report(' Rotation has not converged !') ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR,A,natu,XS,XB,RMS) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR,A,natu,XS,XB,RMS) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR,A,natu,XS,XB,RMS) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR,A,natu,XS,XB,RMS) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(ANG,A,natu,XS,XB,RMS) implicit none INTEGER*4 NATu,I REAL*8 S1,S2,S3,C1,C2,C3,TX,TY,TZ,RMS,ANG(3),FU real*8 A(3,3),XS(natu,3),XB(natu,3) S1=SIN(ANG(1)) S2=SIN(ANG(2)) S3=SIN(ANG(3)) C1=COS(ANG(1)) C2=COS(ANG(2)) C3=COS(ANG(3)) C ang1 = theta C ang2 = phi C ang3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 RMS=0.0d0 DO 1 I=1,NATu TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 RMS=RMS+(TX-XB(I,1))*(TX-XB(I,1))+ 1 (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3)) FU=RMS RETURN END SUBROUTINE READTEN(nat,P,A,LAPT,x) IMPLICIT none integer*4 nat,L,I,J,ID,IG REAL*8 P(nat,3,3),BOHR,A(nat,3,3),x(*),SUM,EPS real*8,allocatable:: R(:,:) LOGICAL LAPT BOHR=0.52917705993d0 allocate(R(3,nat)) DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I+3*(L-1))/BOHR open(15,file='FILE.TEN',status='old') READ (15,*) C C Read dipole derivatives: DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) C C If no axial tensor, then construct its polar part only: IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SUM=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2203 A(L,J,I)=SUM else DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) endif WRITE(6,*)' VCD and IR parameters read in ...' C C Put axial tensor in the origin on the atom L: DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(6,*)' AAT shifted to atomic origin' close(15) RETURN END SUBROUTINE tenback(nat,P,A,x) IMPLICIT none integer*4 nat,L,I,J,ID,IG REAL*8 P(nat,3,3),BOHR,A(nat,3,3),x(*),SUM,EPS,R(3) BOHR=0.52917705993d0 C Put axial tensor in common origin: DO 2201 L=1,NAT R(1)=X(1+3*(L-1))/BOHR R(2)=X(2+3*(L-1))/BOHR R(3)=X(3+3*(L-1))/BOHR DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(6,*)' AAT shifted to common origin' RETURN END FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 RETURN END subroutine readopt(uu,nat,z,nu,natu,m1,m2,N,iwr,LAPT,qmax, 1lsmat,wmin,wmax,delq,umax,umin,qvmi,qvma,ltab,lprn, 1np,wsmin,wsmax,fwhm,lrot,lqpt,ldip,lg,lx,lttt,ENM,temp,cpol, 1nd,gau,vmin,vmax,v,ux,uy,epsr,is,tdc) IMPLICIT none character*4 S4 integer*4 nu,nat,natu,uu(nat,nat),N,m1,m2,z(*),i,j,ia,iwr, 1umax,umin,np,nd,is real*8 wmin,wmax,delq,qvma,qvmi,qmax,wsmin,wsmax,fwhm,ENM,temp, 1gau,gamma,vmax,vmin,vcm(4),v(4),ux(*),uy(*),epsr,CM logical LAPT,lsmat,ltab,lprn,lrot,lqpt,ldip,lg,lx,lttt,cpol, 1ldeut,tdc CM=219474.63d0 nu=0 natu=0 N=0 m1=0 m2=0 iwr=0 c maximal q-vector (if zero then set N-1) qmax=0.0d0 c correction for vibrationalpolarizability cpol=.false. c involved units in cpol -nd ... nd nd=3 c bandwidth for vibr polarizability gamma=20.0d0 c frequency limits for teh vibr p[ol correction vmin=1400.0d0 vmax=1800.0d0 c vibrational permittivity epsr=3.0d0 c transition dipole compling correction tdc=.false. c deuterated amide: ldeut=.false. c Amide I H: c vibrational frequency in cm-1: vcm(1)=1675.0d0 c projection of dipole to x ~ C->N bond direction c z ~ CN x CO, y ~ z x x ux(1)=-0.101d0 c projection of dipole to y-direction: uy(1)=+0.103d0 c Amide II H: vcm(2)=1552.0d0 ux(2)=-0.117d0 uy(2)=-0.024d0 c Amide I D: vcm(3)=1667.0d0 ux(3)=-0.116d0 uy(3)=+0.100d0 c Amide II D: vcm(4)=1482.0d0 ux(4)= 0.058d0 uy(4)=+0.054d0 LAPT=.false. c print S-matrix lsmat=.false. c rotation using complex coordinates: lrot=.false. c read quadrupole derivatives from FILE.QEN lqpt=.false. c write dipoles for all transitions to DIP.TXT: ldip=.false. c write control coordinates to BIGC.X: lx=.false. c Gaussian bands: lg=.false. c Do also Raman and ROA: lttt=.true. c Raman excitation wavelength: ENM=532.0d0 c temperature in K: temp=293.0d0 c minimal and maximal q-vector for lsmat qvmi=0.0d0 qvma=0.0d0 c minimal and maximal wavenumber for lsmat: wmin=0.0d0 wmax=0.0d0 delq=1.0d0 c print tables with all transitions: ltab=.true. c print smooth spectra: lprn=.true. c number of points in the spectra: np=3951 wsmin=50.0d0 wsmax=4000.0d0 c FUll width at half maximum: fwhm=10.0d0 c maximal and minimal unit for LSMAT umax=0 umin=0 N=0 c NUNI: c nu natu number of input units and number of atoms in each c uu the atomic list for each unit c MAIN c m1 m2 the amin (reference) unit and a neighboring one c LENGTH length/number of units in the final structure c QMAX maximum wave vector c DELQ increment of the wave vector c WRITE controls amount of output c APT AAT from APT open(2,file='CPROPAG.OPT',status='old') 1 read(2,80,end=88,err=88)s4 80 format(a4) if(s4.eq.'NUNI')then read(2,*)nu,natu if(nu.gt.nat)call report('nu>nat!') if(natu.gt.nat)call report('natu>nat!') do 2 i=1,nu read(2,*)(uu(i,ia),ia=1,natu) do 2 j=1,i-1 do 21 ia=1,natu if( z(uu(i,ia)).ne.z(uu(j,ia)) )then write(6,*)' units ',i,j,', atom ',ia,z(uu(i,ia)),z(uu(j,ia)) call report('different types') endif 21 continue 2 continue endif if(s4.eq.'MAIN')read(2,*)m1,m2 if(s4.eq.'LENG')read(2,*)N if(s4.eq.'DELQ')read(2,*)delq if(s4.eq.'DEUT')read(2,*)ldeut if(s4.eq.'EPSR')read(2,*)epsr if(s4.eq.'FWHM')read(2,*)fwhm if(s4.eq.'LAMB')read(2,*)ENM if(s4.eq.'LAPT')read(2,*)LAPT if(s4.eq.'CPOL')read(2,*)cpol if(s4.eq.'LTDC')read(2,*)tdc if(s4.eq.'GAMM')read(2,*)gamma if(s4.eq.'LDIP')read(2,*)ldip if(s4.eq.'LGAU')read(2,*)lg if(s4.eq.'LPRN')read(2,*)lprn if(s4.eq.'LQPT')read(2,*)lqpt if(s4.eq.'LROT')read(2,*)lrot if(s4.eq.'LSMA')read(2,*)lsmat if(s4.eq.'LTAB')read(2,*)ltab if(s4.eq.'LTTT')read(2,*)lttt if(s4.eq.'LXCO')read(2,*)lx if(s4.eq.'NDND')read(2,*)nd if(s4.eq.'NPOI')read(2,*)np if(s4.eq.'QMAX')read(2,*)qmax if(s4.eq.'QVMA')read(2,*)qvma if(s4.eq.'QVMI')read(2,*)qvmi if(s4.eq.'TEMP')read(2,*)temp if(s4.eq.'UMAX')read(2,*)umax if(s4.eq.'UMIN')read(2,*)umin if(s4.eq.'WMAX')read(2,*)wmax if(s4.eq.'WMIN')read(2,*)wmin if(s4.eq.'WRIT')read(2,*)iwr if(s4.eq.'WSMA')read(2,*)wsmax if(s4.eq.'WSMI')read(2,*)wsmin if(s4.eq.'VMAX')read(2,*)vmax if(s4.eq.'VMIN')read(2,*)vmin goto 1 88 close(2) if(N.eq.0)call report('N not set') if(nu.eq.0)call report('NUNI not set') if(natu.eq.0)call report('natu not set') if(natu.lt.3)call report('natu < 3') if(m1.eq.0)call report('MAIN 1 not set') if(m2.eq.0)call report('MAIN 2 not set') if(dabs(qmax).lt.1.0d-10)qmax=dble(N-1) if(lrot.and.lsmat)call report('S-matrix wiht LROT not possible') umin=max(umin,1) umax=min(umax,N) gau=gamma/CM if(ldeut)then is=2 else is=0 endif do 3 i=1,4 3 v(i)=vcm(i)/CM return end SUBROUTINE report(s) character*(*) s write(6,*)s stop end SUBROUTINE READX(ic,nat,z,r) IMPLICIT none real*8 r(*) integer*4 ic,nat,i,ix,z(*) open(2,file='FILE.X',status='old') READ(2,*) read(2,*)nat if(ic.eq.1)then do 1 i=1,nat 1 READ(2,*)z(i),(r(ix+3*(i-1)),ix=1,3) endif close(2) RETURN END c ============================================================ subroutine tql2(nm,n,d,e,z,ierr) c integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr double precision d(n),e(n),z(nm,n) double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag c c this subroutine is a translation of the algol procedure tql2, c num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and c wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). c c this subroutine finds the eigenvalues and eigenvectors c of a symmetric tridiagonal matrix by the ql method. c the eigenvectors of a full symmetric matrix can also c be found if tred2 has been used to reduce this c full matrix to tridiagonal form. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e contains the subdiagonal elements of the input matrix c in its last n-1 positions. e(1) is arbitrary. c c z contains the transformation matrix produced in the c reduction by tred2, if performed. if the eigenvectors c of the tridiagonal matrix are desired, z must contain c the identity matrix. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct but c unordered for indices 1,2,...,ierr-1. c c e has been destroyed. c c z contains orthonormal eigenvectors of the symmetric c tridiagonal (or full) matrix. if an error exit is made, c z contains the eigenvectors associated with the stored c eigenvalues. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e(i-1) = e(i) c f = 0.0d0 tst1 = 0.0d0 e(n) = 0.0d0 c do 240 l = 1, n j = 0 h = dabs(d(l)) + dabs(e(l)) if (tst1 .lt. h) tst1 = h c .......... look for small sub-diagonal element .......... do 110 m = l, n tst2 = tst1 + dabs(e(m)) if (tst2 .eq. tst1) go to 120 c .......... e(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 220 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 l2 = l1 + 1 g = d(l) p = (d(l1) - g) / (2.0d0 * e(l)) r = pythag(p,1.0d0) d(l) = e(l) / (p + dsign(r,p)) d(l1) = e(l) * (p + dsign(r,p)) dl1 = d(l1) h = g - d(l) if (l2 .gt. n) go to 145 c do 140 i = l2, n 140 d(i) = d(i) - h c 145 f = f + h c .......... ql transformation .......... p = d(m) c = 1.0d0 c2 = c el1 = e(l1) s = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml c3 = c2 c2 = c s2 = s i = m - ii g = c * e(i) h = c * p r = pythag(p,e(i)) e(i+1) = s * r s = e(i) / r c = p / r p = c * d(i) - s * g d(i+1) = h + s * (c * g + s * d(i)) c .......... form vector .......... do 180 k = 1, n h = z(k,i+1) if(dabs(z(k,i)).lt.1.0d-40)z(k,i)=0.0d0 z(k,i+1) = s * z(k,i) + c * h z(k,i) = c * z(k,i) - s * h 180 continue c 200 continue c p = -s * s2 * c3 * el1 * e(l) / dl1 e(l) = s * p d(l) = c * p tst2 = tst1 + dabs(e(l)) if (tst2 .gt. tst1) go to 130 220 d(l) = d(l) + f 240 continue c .......... order eigenvalues and eigenvectors .......... do 300 ii = 2, n i = ii - 1 k = i p = d(i) c do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue c if (k .eq. i) go to 300 d(k) = d(i) d(i) = p c do 280 j = 1, n p = z(j,i) z(j,i) = z(j,k) z(j,k) = p 280 continue c 300 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end c ============================================================ subroutine htridi(nm,n,ar,ai,d,e,e2,tau) c integer i,j,k,l,n,ii,nm,jp1 double precision ar(nm,n),ai(nm,n),d(n),e(n),e2(n),tau(2,n) double precision f,g,h,fi,gi,hh,si,scale,pythag c c this subroutine is a translation of a complex analogue of c the algol procedure tred1, num. math. 11, 181-195(1968) c by martin, reinsch, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). c c this subroutine reduces a complex hermitian matrix c to a real symmetric tridiagonal matrix using c unitary similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c ar and ai contain the real and imaginary parts, c respectively, of the complex hermitian input matrix. c only the lower triangle of the matrix need be supplied. c c on output c c ar and ai contain information about the unitary trans- c formations used in the reduction in their full lower c triangles. their strict upper triangles and the c diagonal of ar are unaltered. c c d contains the diagonal elements of the the tridiagonal matrix. c c e contains the subdiagonal elements of the tridiagonal c matrix in its last n-1 positions. e(1) is set to zero. c c e2 contains the squares of the corresponding elements of e. c e2 may coincide with e if the squares are not needed. c c tau contains further information about the transformations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c tau(1,n) = 1.0d0 tau(2,n) = 0.0d0 c do 100 i = 1, n 100 d(i) = ar(i,i) c .......... for i=n step -1 until 1 do -- .......... do 300 ii = 1, n i = n + 1 - ii l = i - 1 h = 0.0d0 scale = 0.0d0 if (l .lt. 1) go to 130 c .......... scale row (algol tol then not needed) .......... do 120 k = 1, l 120 scale = scale + dabs(ar(i,k)) + dabs(ai(i,k)) c if (scale .ne. 0.0d0) go to 140 tau(1,l) = 1.0d0 tau(2,l) = 0.0d0 130 e(i) = 0.0d0 e2(i) = 0.0d0 go to 290 c 140 do 150 k = 1, l ar(i,k) = ar(i,k) / scale ai(i,k) = ai(i,k) / scale h = h + ar(i,k) * ar(i,k) + ai(i,k) * ai(i,k) 150 continue c e2(i) = scale * scale * h g = dsqrt(h) e(i) = scale * g f = pythag(ar(i,l),ai(i,l)) c .......... form next diagonal element of matrix t .......... if (f .eq. 0.0d0) go to 160 tau(1,l) = (ai(i,l) * tau(2,i) - ar(i,l) * tau(1,i)) / f si = (ar(i,l) * tau(2,i) + ai(i,l) * tau(1,i)) / f h = h + f * g g = 1.0d0 + g / f ar(i,l) = g * ar(i,l) ai(i,l) = g * ai(i,l) if (l .eq. 1) go to 270 go to 170 160 tau(1,l) = -tau(1,i) si = tau(2,i) ar(i,l) = g 170 f = 0.0d0 c do 240 j = 1, l g = 0.0d0 gi = 0.0d0 c .......... form element of a*u .......... do 180 k = 1, j g = g + ar(j,k) * ar(i,k) + ai(j,k) * ai(i,k) gi = gi - ar(j,k) * ai(i,k) + ai(j,k) * ar(i,k) 180 continue c jp1 = j + 1 if (l .lt. jp1) go to 220 c do 200 k = jp1, l g = g + ar(k,j) * ar(i,k) - ai(k,j) * ai(i,k) gi = gi - ar(k,j) * ai(i,k) - ai(k,j) * ar(i,k) 200 continue c .......... form element of p .......... 220 e(j) = g / h tau(2,j) = gi / h f = f + e(j) * ar(i,j) - tau(2,j) * ai(i,j) 240 continue c hh = f / (h + h) c .......... form reduced a .......... do 260 j = 1, l f = ar(i,j) g = e(j) - hh * f e(j) = g fi = -ai(i,j) gi = tau(2,j) - hh * fi tau(2,j) = -gi c do 260 k = 1, j ar(j,k) = ar(j,k) - f * e(k) - g * ar(i,k) x + fi * tau(2,k) + gi * ai(i,k) ai(j,k) = ai(j,k) - f * tau(2,k) - g * ai(i,k) x - fi * e(k) - gi * ar(i,k) 260 continue c 270 do 280 k = 1, l ar(i,k) = scale * ar(i,k) ai(i,k) = scale * ai(i,k) 280 continue c tau(2,l) = -si 290 hh = d(i) d(i) = ar(i,i) ar(i,i) = hh ai(i,i) = scale * dsqrt(h) 300 continue c return end c ============================================================ subroutine tqlrat(n,d,e2,ierr) c integer i,j,l,m,n,ii,l1,mml,ierr double precision d(n),e2(n) double precision b,c,f,g,h,p,r,s,t,epslon,pythag c c this subroutine is a translation of the algol procedure tqlrat, c algorithm 464, comm. acm 16, 689(1973) by reinsch. c c this subroutine finds the eigenvalues of a symmetric c tridiagonal matrix by the rational ql method. c c on input c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e2 contains the squares of the subdiagonal elements of the c input matrix in its last n-1 positions. e2(1) is arbitrary. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct and c ordered for indices 1,2,...ierr-1, but may not be c the smallest eigenvalues. c c e2 has been destroyed. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e2(i-1) = e2(i) c f = 0.0d0 t = 0.0d0 e2(n) = 0.0d0 c do 290 l = 1, n j = 0 h = dabs(d(l)) + dsqrt(e2(l)) if (t .gt. h) go to 105 t = h b = epslon(t) c = b * b c .......... look for small squared sub-diagonal element .......... 105 do 110 m = l, n if (e2(m) .le. c) go to 120 c .......... e2(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 210 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 s = dsqrt(e2(l)) g = d(l) p = (d(l1) - g) / (2.0d0 * s) r = pythag(p,1.0d0) d(l) = s / (p + dsign(r,p)) h = g - d(l) c do 140 i = l1, n 140 d(i) = d(i) - h c f = f + h c .......... rational ql transformation .......... g = d(m) if (g .eq. 0.0d0) g = b h = g s = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml i = m - ii p = g * h r = p + e2(i) e2(i+1) = s * r s = e2(i) / r d(i+1) = h + s * (h + d(i)) g = d(i) - e2(i) / g if (g .eq. 0.0d0) g = b h = g * p / r 200 continue c e2(l) = s * g d(l) = h c .......... guard against underflow in convergence test .......... if (h .eq. 0.0d0) go to 210 if (dabs(e2(l)) .le. dabs(c/h)) go to 210 e2(l) = h * e2(l) if (e2(l) .ne. 0.0d0) go to 130 210 p = d(l) + f c .......... order eigenvalues .......... if (l .eq. 1) go to 250 c .......... for i=l step -1 until 2 do -- .......... do 230 ii = 2, l i = l + 2 - ii if (p .ge. d(i-1)) go to 270 d(i) = d(i-1) 230 continue c 250 i = 1 270 d(i) = p 290 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end c ============================================================ subroutine htribk(nm,n,ar,ai,tau,m,zr,zi) c integer i,j,k,l,m,n,nm double precision ar(nm,n),ai(nm,n),tau(2,n),zr(nm,m),zi(nm,m) double precision h,s,si c c this subroutine is a translation of a complex analogue of c the algol procedure trbak1, num. math. 11, 181-195(1968) c by martin, reinsch, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). c c this subroutine forms the eigenvectors of a complex hermitian c matrix by back transforming those of the corresponding c real symmetric tridiagonal matrix determined by htridi. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c ar and ai contain information about the unitary trans- c formations used in the reduction by htridi in their c full lower triangles except for the diagonal of ar. c c tau contains further information about the transformations. c c m is the number of eigenvectors to be back transformed. c c zr contains the eigenvectors to be back transformed c in its first m columns. c c on output c c zr and zi contain the real and imaginary parts, c respectively, of the transformed eigenvectors c in their first m columns. c c note that the last component of each returned vector c is real and that vector euclidean norms are preserved. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c if (m .eq. 0) go to 200 c .......... transform the eigenvectors of the real symmetric c tridiagonal matrix to those of the hermitian c tridiagonal matrix. .......... do 50 k = 1, n c do 50 j = 1, m zi(k,j) = -zr(k,j) * tau(2,k) zr(k,j) = zr(k,j) * tau(1,k) 50 continue c if (n .eq. 1) go to 200 c .......... recover and apply the householder matrices .......... do 140 i = 2, n l = i - 1 h = ai(i,i) if (h .eq. 0.0d0) go to 140 c do 130 j = 1, m s = 0.0d0 si = 0.0d0 c do 110 k = 1, l s = s + ar(i,k) * zr(k,j) - ai(i,k) * zi(k,j) si = si + ar(i,k) * zi(k,j) + ai(i,k) * zr(k,j) 110 continue c .......... double divisions avoid possible underflow .......... s = (s / h) / h si = (si / h) / h c do 120 k = 1, l zr(k,j) = zr(k,j) - s * ar(i,k) - si * ai(i,k) zi(k,j) = zi(k,j) - si * ar(i,k) + s * ai(i,k) 120 continue c 130 continue c 140 continue c 200 return end c ============================================================ subroutine ch(n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) c integer i,j,n,nm,ierr,matz double precision ar(n,n),ai(n,n),w(n),zr(n,n),zi(n,n), x fv1(n),fv2(n),fm1(2,n) c c this subroutine calls the recommended sequence of c subroutines from the eigensystem subroutine package (eispack) c to find the eigenvalues and eigenvectors (if desired) c of a complex hermitian matrix. c c on input c c Petr Bour Change: nm=0 by default! c c nm must be set to the row dimension of the two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix a=(ar,ai). c c ar and ai contain the real and imaginary parts, c respectively, of the complex hermitian matrix. c c matz is an integer variable set equal to zero if c only eigenvalues are desired. otherwise it is set to c any non-zero integer for both eigenvalues and eigenvectors. c c on output c c w contains the eigenvalues in ascending order. c c zr and zi contain the real and imaginary parts, c respectively, of the eigenvectors if matz is not zero. c c ierr is an integer output variable set equal to an error c completion code described in the documentation for tqlrat c and tql2. the normal completion code is zero. c c fv1, fv2, and fm1 are temporary storage arrays. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c nm=n if (n .le. nm) go to 10 ierr = 10 * n go to 50 c 10 call htridi(nm,n,ar,ai,w,fv1,fv2,fm1) if (matz .ne. 0) go to 20 c .......... find eigenvalues only .......... call tqlrat(n,w,fv2,ierr) go to 50 c .......... find both eigenvalues and eigenvectors .......... 20 do 40 i = 1, n c do 30 j = 1, n zr(j,i) = 0.0d0 30 continue c zr(i,i) = 1.0d0 40 continue c call tql2(nm,n,w,fv1,zr,ierr) if (ierr .ne. 0) go to 50 call htribk(nm,n,ar,ai,fm1,n,zr,zi) 50 return end double precision function epslon (x) double precision x c c estimate unit roundoff in quantities of size x. c double precision a,b,c,eps c c this program should function properly on all systems c satisfying the following two assumptions, c 1. the base used in representing floating point c numbers is not a power of three. c 2. the quantity a in statement 10 is represented to c the accuracy used in floating point variables c that are stored in memory. c the statement number 10 and the go to 10 are intended to c force optimizing compilers to generate code satisfying c assumption 2. c under these assumptions, it should be true that, c a is not exactly equal to four-thirds, c b has a zero for its last bit or digit, c c is not exactly equal to one, c eps measures the separation of 1.0 from c the next larger floating point number. c the developers of eispack would appreciate being informed c about any systems where these assumptions do not hold. c c this version dated 4/6/83. c a = 4.0d0/3.0d0 10 b = a - 1.0d0 c = b + b + b eps = dabs(c-1.0d0) if (eps .eq. 0.0d0) go to 10 epslon = eps*dabs(x) return end c ============================================================ double precision function pythag(a,b) double precision a,b c c finds dsqrt(a**2+b**2) without overflow or destructive underflow c double precision p,r,s,t,u p = dmax1(dabs(a),dabs(b)) if (p .eq. 0.0d0) go to 20 r = (dmin1(dabs(a),dabs(b))/p)**2 10 continue t = 4.0d0 + r if (t .eq. 4.0d0) go to 20 s = r/t u = 1.0d0 + 2.0d0*s p = u*p r = (s/u)**2 * r go to 10 20 pythag = p return end subroutine pa1(m1,nat,natu,z,u,r,T1,zz,xt) c put atoms of ref. unit to xt,zz, relative to center T1 implicit none integer*4 m1,nat,natu,z(*),zz(natu),ia,ix,u(nat,nat) real*8 r(*),xt(3*natu),T1(3) do 3 ia=1,natu zz(ia)=z(u(m1,ia)) do 3 ix=1,3 3 xt(ix+3*(ia-1))=r(ix+3*(u(m1,ia)-1))-T1(ix) return end c ===========================================================C subroutine gengeo(N,nat,nu3,natu,v,A,r,m1,u,z,T1) implicit none real*8 v(*),A(3,3),r(*),T1(3),rr(3) integer*4 nat,m1,N,u(nat,nat),natu,nu3,z(*),ia,I,ix,jx,ii real*8,allocatable::xt(:),rrf(:) integer*4,allocatable::zz(:) allocate(xt(nu3),zz(natu),rrf(nu3)) call pa1(m1,nat,natu,z,u,r,T1,zz,xt) open(9,file='BIGC.X') write(9,*)'1D crystal geometry' write(9,*)natu*N rrf=xt c initial geometry: do 4 ia=1,natu 4 write(9,600)zz(ia),(T1(ix)+xt(ix+3*(ia-1)),ix=1,3) 600 format(i4,3f12.6,' 0 0 0 0 0 0 0 0.0') c loop over other cells: do 1 I=1,N-1 c shift cell forward c rotate previous geometry in rrf: xt=rrf do 5 ia=1,natu do 5 ix=1,3 ii=ix+3*(ia-1) rrf(ii)=0.0d0 do 5 jx=1,3 5 rrf(ii)=rrf(ii)+A(ix,jx)*xt(jx+3*(ia-1)) do 51 ia=1,natu do 52 ix=1,3 52 rr(ix)=T1(ix)+dble(I)*v(ix)+rrf(ix+3*(ia-1)) 51 write(9,600)zz(ia),(rr(ix),ix=1,3) 1 continue close(9) write(6,*)'BIGC.X written' return end subroutine trs(natu,u,A) implicit none integer*4 natu,i real*8 u(*),A(3,3),x,y,z do 1 i=1,natu x=u(1+3*(i-1)) y=u(2+3*(i-1)) z=u(3+3*(i-1)) u(1+3*(i-1))=A(1,1)*x+A(1,2)*y+A(1,3)*z u(2+3*(i-1))=A(2,1)*x+A(2,2)*y+A(2,3)*z 1 u(3+3*(i-1))=A(3,1)*x+A(3,2)*y+A(3,3)*z return end subroutine putaat(Xcell,Pat,Aat,Acom,natu) implicit none integer*4 natu,L,J,I,ID,IG,ii real*8 Xcell(*),Pat(natu*3,3),Aat(natu*3,3),Acom(natu*3,3),SUM, 1EPS DO 2201 L=1,natu DO 2201 J=1,3 ii=J+3*(L-1) DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM+0.25d0*EPS(I,IG,ID)*Xcell(IG+3*(L-1))*Pat(ii,ID) 2201 Acom(ii,I)=Aat(ii,I)+SUM return end subroutine putqat(Xcell,P,Qat,Qcom,natu) implicit none integer*4 natu,L,I,LI real*8 Xcell(*),P(natu*3,3),Qat(natu*3,6),Qcom(natu*3,6),R(3) DO 1 L=1,NATu do 2 I=1,3 2 R(I)=Xcell(I+3*(L-1)) DO 1 I=1,3 LI=I+3*(L-1) Qcom(LI,1)=Qat(LI,1)+(P(LI,1)*R(1)+P(LI,1)*R(1)) Qcom(LI,2)=Qat(LI,2)+(P(LI,2)*R(2)+P(LI,2)*R(2)) Qcom(LI,3)=Qat(LI,3)+(P(LI,3)*R(3)+P(LI,3)*R(3)) Qcom(LI,4)=Qat(LI,4)+(P(LI,1)*R(2)+P(LI,2)*R(1)) Qcom(LI,5)=Qat(LI,5)+(P(LI,1)*R(3)+P(LI,3)*R(1)) 1 Qcom(LI,6)=Qat(LI,6)+(P(LI,2)*R(3)+P(LI,3)*R(2)) return end subroutine trr(Pat,natu,am) implicit none integer*4 natu,i,ix,jx,ii,jj real*8 Pat(natu*3,3),t(3,3),am(3,3) do 1 i=1,natu do 11 ix=1,3 do 11 jj=1,3 t(ix,jj)=0.0d0 do 11 jx=1,3 11 t(ix,jj)=t(ix,jj)+am(jj,jx)*Pat(3*(i-1)+ix,jx) do 1 ii=1,3 do 1 jj=1,3 Pat(3*(i-1)+ii,jj)=0.0d0 do 1 ix=1,3 1 Pat(3*(i-1)+ii,jj)=Pat(3*(i-1)+ii,jj)+am(ii,ix)*t(ix,jj) return end subroutine tra(a,natu,am) implicit none integer*4 natu,i,j,k,l real*8 a(natu*3,3,3),t(3,3,3),u(3,3,3),am(3,3) do 2 i=1,natu do 1 j=1,3 do 1 k=1,3 do 1 l=1,3 1 t(j,k,l)=am(j,1)*a(3*(i-1)+1,k,l) 1 +am(j,2)*a(3*(i-1)+2,k,l) 1 +am(j,3)*a(3*(i-1)+3,k,l) do 3 j=1,3 do 3 k=1,3 do 3 l=1,3 3 u(j,k,l)=+am(k,1)*t(j,1,l)+am(k,2)*t(j,2,l)+am(k,3)*t(j,3,l) do 2 j=1,3 do 2 k=1,3 do 2 l=1,3 2 a(3*(i-1)+j,k,l)=am(l,1)*u(j,k,1) 1 +am(l,2)*u(j,k,2) 1 +am(l,3)*u(j,k,3) return end subroutine trt(a,natu,am) implicit none integer*4 natu,i,jx,kx,lx,ii,jj,kk,ll real*8 a(natu*3,3,3,3),t(3,3,3,3),u(3,3,3,3),am(3,3) do 2 i=1,natu do 1 ii=1,3 do 1 jx=1,3 do 1 kx=1,3 do 1 lx=1,3 1 t(ii,jx,kx,lx)=am(ii,1)*a(3*(i-1)+1,jx,kx,lx) 1 +am(ii,2)*a(3*(i-1)+2,jx,kx,lx) 1 +am(ii,3)*a(3*(i-1)+3,jx,kx,lx) do 3 ii=1,3 do 3 jj=1,3 do 3 kx=1,3 do 3 lx=1,3 3 u(ii,jj,kx,lx)=am(jj,1)*t(ii,1,kx,lx) 1 +am(jj,2)*t(ii,2,kx,lx) 1 +am(jj,3)*t(ii,3,kx,lx) do 4 ii=1,3 do 4 jj=1,3 do 4 kk=1,3 do 4 lx=1,3 4 t(ii,jj,kk,lx)=am(kk,1)*u(ii,jj,1,lx) 1 +am(kk,2)*u(ii,jj,2,lx) 1 +am(kk,3)*u(ii,jj,3,lx) do 2 ii=1,3 do 2 jj=1,3 do 2 kk=1,3 do 2 ll=1,3 2 a(3*(i-1)+ii,jj,kk,ll)=am(ll,1)*t(ii,jj,kk,1) 1 +am(ll,2)*t(ii,jj,kk,2) 1 +am(ll,3)*t(ii,jj,kk,3) return end subroutine trq(Qat,natu,am) implicit none integer*4 natu,i,ix,jx,ii,jj,kk,kx real*8 Qat(natu*3,6),t(3,3,3),am(3,3),Q(3,3,3) do 1 i=1,natu do 11 ix=1,3 Q(ix,1,1)=Qat(3*(i-1)+ix,1) Q(ix,2,2)=Qat(3*(i-1)+ix,2) Q(ix,3,3)=Qat(3*(i-1)+ix,3) Q(ix,1,2)=Qat(3*(i-1)+ix,4) Q(ix,2,1)=Qat(3*(i-1)+ix,4) Q(ix,1,3)=Qat(3*(i-1)+ix,5) Q(ix,3,1)=Qat(3*(i-1)+ix,5) Q(ix,2,3)=Qat(3*(i-1)+ix,6) 11 Q(ix,3,2)=Qat(3*(i-1)+ix,6) do 2 ix=1,3 do 2 jx=1,3 do 2 kk=1,3 t(ix,jx,kk)=0.0d0 do 2 kx=1,3 2 t(ix,jx,kk)=t(ix,jx,kk)+am(kk,kx)*Q(ix,jx,kx) do 3 ix=1,3 do 3 jj=1,3 do 3 kk=1,3 Q(ix,jj,kk)=0.0d0 do 3 jx=1,3 3 Q(ix,jj,kk)=Q(ix,jj,kk)+am(jj,jx)*t(ix,jx,kk) do 4 ii=1,3 do 4 jj=1,3 do 4 kk=1,3 t(ii,jj,kk)=0.0d0 do 4 ix=1,3 4 t(ii,jj,kk)=t(ii,jj,kk)+am(ii,ix)*Q(ix,jj,kk) do 1 ix=1,3 Qat(3*(i-1)+ix,1)=t(ix,1,1) Qat(3*(i-1)+ix,2)=t(ix,2,2) Qat(3*(i-1)+ix,3)=t(ix,3,3) Qat(3*(i-1)+ix,4)=t(ix,1,2) Qat(3*(i-1)+ix,5)=t(ix,1,3) 1 Qat(3*(i-1)+ix,6)=t(ix,3,2) return end subroutine sm(sr,si,m,n) implicit none integer*4 n,i,j,ia real*8 sr(n,n),si(n,n),m(*) do 1 i=1,n ia=(i+2)/3 do 1 j=1,n si(i,j)=si(i,j)/dsqrt(m(ia)) 1 sr(i,j)=sr(i,j)/dsqrt(m(ia)) return end subroutine trc(Xcell,am,natu,v) c Xcell in bohrs and v in A implicit none real*8 Xcell(*),am(3,3),mc(3),v(*),BOHR integer*4 ix,ia,natu real*8,allocatable::rt(:) allocate(rt(3*natu)) BOHR=0.52917705993d0 c mass center: do 18 ix=1,3 mc(ix)=0.0d0 do 151 ia=1,natu 151 mc(ix)=mc(ix)+Xcell(ix+3*(ia-1)) 18 mc(ix)=mc(ix)/dble(natu) c subtract mass center: do 16 ia=1,natu do 16 ix=1,3 16 rt(ix+3*(ia-1))=Xcell(ix+3*(ia-1))-mc(ix) c rotate relative coordinates with supplied matrix: call trs(natu,rt,am) c shift and add to mass center: do 17 ia=1,natu do 17 ix=1,3 17 Xcell(ix+3*(ia-1))=mc(ix)+v(ix)/BOHR+rt(ix+3*(ia-1)) return end subroutine wrs(natu,lsmat,W,rq,wmin,wmax,qvmi,qvma, 1itt,it,I,svr,svi,umin,umax,sI,cI) implicit none logical lsmat integer*4 itt,it,ia,ix,I,natu,umin,umax,jj real*8 W,rq,wmin,wmax,qvmi,qvma,svr(*),svi(*),prd, 1sI,cI,re,im prd=1.0d-4 if(lsmat.and.W.ge.wmin.and.W.le.wmax.and. 1 I+1.ge.umin.and.I+1.le.umax. 1 and.rq.ge.qvmi-prd.and.rq.le.qvma+prd)then itt=itt+1 write(46,403)itt,I,it,W 403 format(3i6,f12.4) do 20 ia=1,natu do 20 ix=1,3 jj=ix+3*(ia-1) re=svr(jj)*cI-svi(jj)*sI im=svr(jj)*sI+svi(jj)*cI 20 write(46,404)ia,ix,re,im 404 format(2i6,3f12.6) endif return end subroutine massa(nat,natu,uu,m,mu,m1) implicit none integer*4 natu,nat,ia,at,uu(nat,nat),m1 real*8 m(*),mu(*) do 1 ia=1,natu at=uu(m1,ia) 1 mu(ia)=m(at) return end c ================================================================== function sf(d,g,x,x0,ic) implicit none logical g integer*4 ic real*8 pi,spi,sf,d,x,x0,dd,ab pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) dd=((x-x0)/d)**2 c write(6,*)d,x,x0,g if(ic.eq.30)then c ORD, real parts, only Lorentz possible: ab=(x-x0)*(x+x0)*x**2*x0**2 sf=ab/(ab**2+(d**2)**2) else c dispersion spectra, imaginary part if(g)then if(dd.lt.20.0d0)then sf=exp(-dd)/d/spi else sf=0.0d0 endif else if(dd.lt.1.0d10)then sf=1.0d0/d/(dd+1.0d0)/pi else sf=0.0d0 endif endif endif return end c ================================================================== subroutine rot1(f,a,n3) implicit none integer*4 n3,nat,ia,i,ja,j,ix,jx,ii,jj real*8 f(n3,n3),a(3,3),t(3,3) nat=n3/3 do 1 ia=1,nat do 1 ja=1,nat c first index, do 11 i=1,3 do 11 jx=1,3 jj=jx+3*(ja-1) t(i,jx)=0.0d0 do 11 ix=1,3 ii=ix+3*(ia-1) 11 t(i,jx)=t(i,jx)+a(ix,i)*f(ii,jj) c second index, do 1 i=1,3 ii=i+3*(ia-1) do 1 j=1,3 jj=j+3*(ja-1) f(ii,jj)=0.0d0 do 1 jx=1,3 1 f(ii,jj)=f(ii,jj)+a(jx,j)*t(i,jx) write(6,*)'Force field rotated' return end subroutine trafo3(q,a) implicit none integer*4 i,j,k real*8 q(3,3,3),t(3,3,3),a(3,3) c first index, do 11 i=1,3 do 11 j=1,3 do 11 k=1,3 11 t(i,j,k)=a(1,i)*q(1,j,k)+a(2,i)*q(2,j,k)+a(3,i)*q(3,j,k) c second index, do 2 i=1,3 do 2 j=1,3 do 2 k=1,3 2 q(i,j,k)=a(1,j)*t(i,1,k)+a(2,j)*t(i,2,k)+a(3,j)*t(i,3,k) c third index, do 3 i=1,3 do 3 j=1,3 do 3 k=1,3 3 t(i,j,k)=a(1,k)*q(i,j,1)+a(2,k)*q(i,j,2)+a(3,k)*q(i,j,3) q=t return end subroutine trafo4(q,a) implicit none integer*4 i,j,k,l real*8 q(3,3,3,3),t(3,3,3,3),a(3,3) do 11 i=1,3 do 11 j=1,3 do 11 k=1,3 do 11 l=1,3 11 t(i,j,k,l)=a(1,i)*q(1,j,k,l)+a(2,i)*q(2,j,k,l)+a(3,i)*q(3,j,k,l) do 2 i=1,3 do 2 j=1,3 do 2 k=1,3 do 2 l=1,3 2 q(i,j,k,l)=a(1,j)*t(i,1,k,l)+a(2,j)*t(i,2,k,l)+a(3,j)*t(i,3,k,l) do 3 i=1,3 do 3 j=1,3 do 3 k=1,3 do 3 l=1,3 3 t(i,j,k,l)=a(1,k)*q(i,j,1,l)+a(2,k)*q(i,j,2,l)+a(3,k)*q(i,j,3,l) do 4 i=1,3 do 4 j=1,3 do 4 k=1,3 do 4 l=1,3 4 q(i,j,k,l)=a(1,l)*t(i,j,k,1)+a(2,l)*t(i,j,k,2)+a(3,l)*t(i,j,k,3) return end subroutine rot3(p,a,n3) implicit none integer*4 n3,nat,ia,ix real*8 p(n3,6),a(3,3),q(3,3,3) nat=n3/3 do 1 ia=1,nat do 111 ix=1,3 q(ix,1,1)=p(ix+3*(ia-1),1) q(ix,2,2)=p(ix+3*(ia-1),2) q(ix,3,3)=p(ix+3*(ia-1),3) q(ix,1,2)=p(ix+3*(ia-1),4) q(ix,2,1)=p(ix+3*(ia-1),4) q(ix,1,3)=p(ix+3*(ia-1),5) q(ix,3,1)=p(ix+3*(ia-1),5) q(ix,2,3)=p(ix+3*(ia-1),6) 111 q(ix,3,2)=p(ix+3*(ia-1),6) call trafo3(q,a) do 1 ix=1,3 p(ix+3*(ia-1),1)=q(ix,1,1) p(ix+3*(ia-1),2)=q(ix,2,2) p(ix+3*(ia-1),3)=q(ix,3,3) p(ix+3*(ia-1),4)=q(ix,1,2) p(ix+3*(ia-1),5)=q(ix,1,3) 1 p(ix+3*(ia-1),6)=q(ix,3,2) write(6,*)'quadrupolar tensor rotated' return end subroutine rot4(p,a,n3) implicit none integer*4 n3,nat,ia,ix,jx,kx real*8 p(n3,3,3),a(3,3),q(3,3,3) nat=n3/3 do 1 ia=1,nat do 111 ix=1,3 do 111 jx=1,3 do 111 kx=1,3 111 q(ix,jx,kx)=p(ix+3*(ia-1),jx,kx) call trafo3(q,a) do 1 ix=1,3 do 1 jx=1,3 do 1 kx=1,3 1 p(ix+3*(ia-1),jx,kx)=q(ix,jx,kx) write(6,*)'polarizability rotated' return end subroutine rot5(p,a,n3) implicit none integer*4 n3,nat,ia,ix,jx,kx,lx real*8 p(n3,3,3,3),a(3,3),q(3,3,3,3) nat=n3/3 do 1 ia=1,nat do 111 ix=1,3 do 111 jx=1,3 do 111 kx=1,3 do 111 lx=1,3 111 q(ix,jx,kx,lx)=p(ix+3*(ia-1),jx,kx,lx) call trafo4(q,a) do 1 ix=1,3 do 1 jx=1,3 do 1 kx=1,3 do 1 lx=1,3 1 p(ix+3*(ia-1),jx,kx,lx)=q(ix,jx,kx,lx) write(6,*)'AA polarizability rotated' return end subroutine rot2(p,a,n3) implicit none integer*4 n3,nat,ia,i,j,ix,jx real*8 p(n3/3,3,3),a(3,3),t(3,3) nat=n3/3 do 1 ia=1,nat t=0.0d0 c first index, do 11 i=1,3 do 11 jx=1,3 do 11 ix=1,3 11 t(i,jx)=t(i,jx)+a(ix,i)*p(ia,ix,jx) c second index, do 1 i=1,3 do 1 j=1,3 p(ia,i,j)=0.0d0 do 1 jx=1,3 1 p(ia,i,j)=p(ia,i,j)+a(jx,j)*t(i,jx) write(6,*)'tensor rotated' return end subroutine rott(f,fr,fi,nat,ur,ui,natu,uu,m1) c transforms APT/ATT in complex coocrdinates c dipole (last) index only implicit none integer*4 nat,at,ia,i,j,jx,natu,uu(nat,nat),m1,ii real*8 f(nat,3,3),fr(natu*3,3),fi(natu*3,3),ur(3,3),ui(3,3), 1t(3,3) do 1 ia=1,natu at=uu(m1,ia) c first index just conserve do 11 i=1,3 do 11 j=1,3 11 t(i,j)=f(at,i,j) c second index transform F . U do 1 i=1,3 ii=i+3*(ia-1) do 1 j=1,3 fr(ii,j)=0.0d0 fi(ii,j)=0.0d0 do 1 jx=1,3 fr(ii,j)=fr(ii,j)+ur(j,jx)*t(i,jx) 1 fi(ii,j)=fi(ii,j)+ui(j,jx)*t(i,jx) write(6,*)'Complex tensor' return end subroutine trafo2c(t,qr,qi,ur,ui) c Q = UU T implicit none integer*4 j,k real*8 t(3,3),qr(3,3),qi(3,3),ur(3,3),ui(3,3), 1tr(3,3),ti(3,3) do 1 j=1,3 do 1 k=1,3 qr(j,k)=ur(j,1)*t(1,k)+ur(j,2)*t(2,k)+ur(j,3)*t(3,k) 1 qi(j,k)=ui(j,1)*t(1,k)+ui(j,2)*t(2,k)+ui(j,3)*t(3,k) do 2 j=1,3 do 2 k=1,3 tr(j,k)=ur(k,1)*qr(j,1)+ur(k,2)*qr(j,2)+ur(k,3)*qr(j,3) 1 -ui(k,1)*qi(j,1)-ui(k,2)*qi(j,2)-ui(k,3)*qi(j,3) 2 ti(j,k)=ui(k,1)*qr(j,1)+ui(k,2)*qr(j,2)+ui(k,3)*qr(j,3) 1 +ur(k,1)*qi(j,1)+ur(k,2)*qi(j,2)+ur(k,3)*qi(j,3) qr=tr qi=ti return end subroutine trafo3c(t,qr,qi,ur,ui) c Q = UUU T implicit none integer*4 j,k,l real*8 t(3,3,3),qr(3,3,3),qi(3,3,3),ur(3,3),ui(3,3), 1tr(3,3,3),ti(3,3,3) do 1 j=1,3 do 1 k=1,3 do 1 l=1,3 qr(j,k,l)=ur(j,1)*t(1,k,l)+ur(j,2)*t(2,k,l)+ur(j,3)*t(3,k,l) 1 qi(j,k,l)=ui(j,1)*t(1,k,l)+ui(j,2)*t(2,k,l)+ui(j,3)*t(3,k,l) do 2 j=1,3 do 2 k=1,3 do 2 l=1,3 tr(j,k,l)=ur(k,1)*qr(j,1,l)+ur(k,2)*qr(j,2,l)+ur(k,3)*qr(j,3,l) 1 -ui(k,1)*qi(j,1,l)-ui(k,2)*qi(j,2,l)-ui(k,3)*qi(j,3,l) 2 ti(j,k,l)=ui(k,1)*qr(j,1,l)+ui(k,2)*qr(j,2,l)+ui(k,3)*qr(j,3,l) 1 +ur(k,1)*qi(j,1,l)+ur(k,2)*qi(j,2,l)+ur(k,3)*qi(j,3,l) do 3 j=1,3 do 3 k=1,3 do 3 l=1,3 qr(j,k,l)=ur(l,1)*tr(j,k,1)+ur(l,2)*tr(j,k,2)+ur(l,3)*tr(j,k,3) 1 -ui(l,1)*ti(j,k,1)-ui(l,2)*ti(j,k,2)-ui(l,3)*ti(j,k,3) 3 qi(j,k,l)=ui(l,1)*tr(j,k,1)+ui(l,2)*tr(j,k,2)+ui(l,3)*tr(j,k,3) 1 +ur(l,1)*ti(j,k,1)+ur(l,2)*ti(j,k,2)+ur(l,3)*ti(j,k,3) return end subroutine rotQ(f,fr,fi,nat,ur,ui,natu,uu,m1) c transforms quadrupole derivatives in complex coocrdinates c derivative index conserve implicit none integer*4 nat,ia,at,j,k,ix,natu,uu(nat,nat),m1,ii real*8 f(3*nat,6),fr(natu*3,3,3),fi(natu*3,3,3),ur(3,3),ui(3,3), 1t(3,3),qr(3,3),qi(3,3) do 1 ia=1,natu at=uu(m1,ia) c first index coord derv conserve do 1 ix=1,3 t(1,1)=f(ix+3*(at-1),1) t(2,2)=f(ix+3*(at-1),2) t(3,3)=f(ix+3*(at-1),3) t(1,2)=f(ix+3*(at-1),4) t(2,1)=f(ix+3*(at-1),4) t(1,3)=f(ix+3*(at-1),5) t(3,1)=f(ix+3*(at-1),5) t(2,3)=f(ix+3*(at-1),6) t(3,2)=f(ix+3*(at-1),6) call trafo2c(t,qr,qi,ur,ui) ii=ix+3*(ia-1) do 1 j=1,3 do 1 k=1,3 fr(ii,j,k)=qr(j,k) 1 fi(ii,j,k)=qi(j,k) write(6,*)'Complex quadrupole' return end subroutine rottt(f,fr,fi,nat,ur,ui,natu,uu,m1) c transforms pol. ders. derivatives in complex coocrdinates c derivative index conserve implicit none integer*4 nat,ia,at,j,k,ix,jx,kx,natu,uu(nat,nat),m1,ii real*8 f(3*nat,3,3),fr(natu*3,3,3),fi(natu*3,3,3),ur(3,3),ui(3,3), 1t(3,3),qr(3,3),qi(3,3) do 1 ia=1,natu at=uu(m1,ia) do 1 ix=1,3 do 2 jx=1,3 do 2 kx=1,3 2 t(jx,kx)=f(ix+3*(at-1),jx,kx) call trafo2c(t,qr,qi,ur,ui) ii=ix+3*(ia-1) do 1 j=1,3 do 1 k=1,3 fr(ii,j,k)=qr(j,k) 1 fi(ii,j,k)=qi(j,k) write(6,*)'Complex polarizability' return end subroutine rott3(f,fr,fi,nat,ur,ui,natu,uu,m1) c transforms pol. ders. derivatives in complex coocrdinates c derivative index conserve implicit none integer*4 nat,ia,at,j,k,ix,jx,kx,natu,uu(nat,nat),m1,ii,l,lx real*8 f(3*nat,3,3,3),fr(natu*3,3,3,3),fi(natu*3,3,3,3), 1ur(3,3),ui(3,3),t(3,3,3),qr(3,3,3),qi(3,3,3) do 1 ia=1,natu at=uu(m1,ia) do 1 ix=1,3 do 2 jx=1,3 do 2 kx=1,3 do 2 lx=1,3 2 t(jx,kx,lx)=f(ix+3*(at-1),jx,kx,lx) call trafo3c(t,qr,qi,ur,ui) ii=ix+3*(ia-1) do 1 j=1,3 do 1 k=1,3 do 1 l=1,3 fr(ii,j,k,l)=qr(j,k,l) 1 fi(ii,j,k,l)=qi(j,k,l) write(6,*)'Complex polarizability' return end subroutine setum(ur,ui) implicit none real*8 ur(3,3),ui(3,3) c set complex transformation matrix c x(complex) = U . x(real) c u - first index new/complex second old/Cartesian c order: x_-1 = (x_x - ix_y)/sqrt(2) c x_0 = x_z c x_+1 = (x_x + ix_y)/sqrt(2) ur=0.0d0 ui=0.0d0 ur(1,1)= 1.0d0/dsqrt(2.0d0) ui(1,2)=-1.0d0/dsqrt(2.0d0) ur(2,3)= 1.0d0 ur(3,1)= 1.0d0/dsqrt(2.0d0) ui(3,2)= 1.0d0/dsqrt(2.0d0) call wm('ur:',ur) call wm('ui:',ui) return end function sp(a,b) real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine norm(v) real*8 v(*),n,sp n=dsqrt(sp(v,v)) v(1)=v(1)/n v(2)=v(2)/n v(3)=v(3)/n return end subroutine tq(qq,nat,x,PP,is) c is=1 : set quadrupole derivatives to local system c is=-1: put them to common origin implicit none integer*4 nat,is,L,I,J,LJ real*8 qq(3*nat,6),x(*),BOHR,PP(nat,3,3),R(3),s BOHR=0.52917705993d0 if(is.eq.1)then s=1.0d0 write(6,*)' QD set to common system' else if(is.eq.-1)then s=-1.0d0 write(6,*)' QD set to local systems' else write(6,*)is call report('tq unknown option') endif endif DO 2204 L=1,nat DO 3 I=1,3 3 R(I)=X(I+3*(L-1))/BOHR DO 2204 J=1,3 LJ=J+3*(L-1) s=1.0d0 qq(LJ,1)=qq(LJ,1)+(PP(L,J,1)*R(1)+PP(L,J,1)*R(1))*s qq(LJ,2)=qq(LJ,2)+(PP(L,J,2)*R(2)+PP(L,J,2)*R(2))*s qq(LJ,3)=qq(LJ,3)+(PP(L,J,3)*R(3)+PP(L,J,3)*R(3))*s qq(LJ,4)=qq(LJ,4)+(PP(L,J,1)*R(2)+PP(L,J,2)*R(1))*s qq(LJ,5)=qq(LJ,5)+(PP(L,J,1)*R(3)+PP(L,J,3)*R(1))*s 2204 qq(LJ,6)=qq(LJ,6)+(PP(L,J,2)*R(3)+PP(L,J,3)*R(2))*s return end subroutine reqd(lqpt,qq,nat,x,PP) implicit none integer*4 i,ia,ix,nat real*8 qq(3*nat,6),x(*),PP(nat,3,3) logical lqpt if(lqpt)then open(9,file='FILE.QEN',status='old') do 1 i=1,9 1 read(9,*) do 2 ia=1,nat do 2 ix=1,3 2 read(9,93)(qq(ix+3*(ia-1),i),i=1,6) 93 format(9x,3e13.4) close(9) write(6,*)'QD read' call tq(qq,nat,x,PP,-1) else qq=0.0d0 write(6,*)'local quadrupole derivatives initiated' endif return end subroutine rewritew(qqr,qr,qqi,qi) real*8 qqr(3,3),qr(6),qqi(3,3),qi(6) qqr(1,1)=qr(1) qqr(2,2)=qr(2) qqr(3,3)=qr(3) qqr(1,2)=qr(4) qqr(2,1)=qr(4) qqr(1,3)=qr(5) qqr(3,1)=qr(5) qqr(2,3)=qr(6) qqr(3,2)=qr(6) qqi(1,1)=qi(1) qqi(2,2)=qi(2) qqi(3,3)=qi(3) qqi(1,2)=qi(4) qqi(2,1)=qi(4) qqi(1,3)=qi(5) qqi(3,1)=qi(5) qqi(2,3)=qi(6) qqi(3,2)=qi(6) return end subroutine trpa(iq,nu3,pr,pi,sr,si,ur,ui) implicit none real*8 sr(nu3,nu3),si(nu3,nu3),ui(3),ur(3), 1pr(nu3,3),pi(nu3,3) integer*4 nu3,iq,k,l do 1 k=1,3 ur(k)=0.0d0 ui(k)=0.0d0 do 1 l=1,nu3 ur(k)=ur(k)+pr(l,k)*sr(l,iq)-pi(l,k)*si(l,iq) 1 ui(k)=ui(k)+pi(l,k)*sr(l,iq)+pr(l,k)*si(l,iq) return end subroutine trpq(iq,nu3,pr,pi,sr,si,ur,ui) implicit none real*8 sr(nu3,nu3),si(nu3,nu3),ui(3,3),ur(3,3), 1pr(nu3,3,3),pi(nu3,3,3) integer*4 nu3,iq,k,l,m do 1 k=1,3 do 1 m=1,3 ur(k,m)=0.0d0 ui(k,m)=0.0d0 do 1 l=1,nu3 ur(k,m)=ur(k,m)+pr(l,k,m)*sr(l,iq)-pi(l,k,m)*si(l,iq) 1 ui(k,m)=ui(k,m)+pi(l,k,m)*sr(l,iq)+pr(l,k,m)*si(l,iq) return end subroutine traa(iq,nu3,pr,pi,sr,si,ur,ui) implicit none real*8 sr(nu3,nu3),si(nu3,nu3),ui(3,3,3),ur(3,3,3), 1pr(nu3,3,3,3),pi(nu3,3,3,3) integer*4 nu3,iq,k,l,m,n do 1 k=1,3 do 1 m=1,3 do 1 n=1,3 ur(k,m,n)=0.0d0 ui(k,m,n)=0.0d0 do 1 l=1,nu3 ur(k,m,n)=ur(k,m,n)+pr(l,k,m,n)*sr(l,iq)-pi(l,k,m,n)*si(l,iq) 1 ui(k,m,n)=ui(k,m,n)+pi(l,k,m,n)*sr(l,iq)+pr(l,k,m,n)*si(l,iq) return end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(*),E(*) REAL*8 MACHEP MACHEP = 1.0d-12 C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END subroutine s1s(k,r,i,N) implicit none real*8 k,r,i,N,down,l l=1.0d0 if(dabs(k).gt.1.0d-5)then c sum(I=0..N-1) 1/sqrt(N) exp[i(2piq/N+kt)I] down=2.0d0*dsqrt(N)*(l-cos(k)) r=(l+cos(k*(N-l))-cos(k)-cos(k*N))/down i=( sin(k*(N-l))+sin(k)-sin(k*N))/down else r=dsqrt(N) i=0.0d0 endif return end subroutine s2s(k,r,i,N) implicit none real*8 k,r,i,N,down,l l=1.0d0 if(dabs(k).gt.1.0d-5)then c sum(I=0..N-1) 1/sqrt(N) I exp[i(2piq/N+kt)I] down=4.0d0*dsqrt(N)*(l-cos(k))**2 r=( (N-l)*cos(k*(N+l))-(3.0d0*N-2.0d0)*cos(k*N) 1 +(3.0d0*N-l)*cos(k*(N-l))-N*cos(k*(N-2.0d0)) 1 +2.0d0*cos(k)-2.0d0)/down i=( (N-l)*sin(k*(N+l))-(3.0d0*N-2.0d0)*sin(k*N) 1 +(3.0d0*N-l)*sin(k*(N-l))-N*sin(k*(N-2.0d0)))/down else r=dsqrt(N)*(N-l)/2.0d0 i=0.0d0 endif return end subroutine srm(a,tau) implicit none real*8 tau,a(3,3) a=0.0d0 a(1,1)= cos(tau) a(1,2)= sin(tau) a(2,1)=-sin(tau) a(2,2)= cos(tau) a(3,3)= 1.0d0 call wm('Rotation matrix in new coordinates',a) return end subroutine rmnc(am,amn,at) implicit none real*8 am(3,3),amn(3,3),at(3,3),t(3,3) integer*4 i,j,ii,jj c AT . AMN: do 1 i=1,3 do 1 j=1,3 t(i,j)=0.0d0 do 1 jj=1,3 1 t(i,j)=t(i,j)+at(i,jj)*amn(jj,j) c AT . AMN . ATt do 2 i=1,3 do 2 j=1,3 am(i,j)=0.0d0 do 2 ii=1,3 2 am(i,j)=am(i,j)+t(i,ii)*at(j,ii) call wm('Rotation matrix in old coordinates',am) return end subroutine wm(s,a) implicit none real*8 a(3,3) integer*4 i,j character*(*)s do 1 i=1,len(s) 1 write(6,601)s(i:i) 601 format(a1,$) write(6,*) write(6,602)((a(i,j),j=1,3),i=1,3) 602 format(3(3f10.4,/)) return end subroutine toQ(mrc,mic,mrt,mit,ar,ai) implicit none integer*4 i,j,k real*8 mrc(3,3),mic(3,3),mrt(3,3),mit(3,3),ar(3,3),ai(3,3), 1ti(3,3),tr(3,3) do 1 i=1,3 do 1 j=1,3 tr(i,j)=0.0d0 ti(i,j)=0.0d0 do 1 k=1,3 c C = U* R tr(i,j)=tr(i,j)+ar(k,i)*mrt(k,j)+ai(k,i)*mit(k,j) 1 ti(i,j)=ti(i,j)-ai(k,i)*mrt(k,j)+ar(k,i)*mit(k,j) do 2 i=1,3 do 2 j=1,3 mrc(i,j)=0.0d0 mic(i,j)=0.0d0 do 2 k=1,3 c C = R U* mrc(i,j)=mrc(i,j)+ar(k,j)*tr(i,k)+ai(k,j)*ti(i,k) 2 mic(i,j)=mic(i,j)-ai(k,j)*tr(i,k)+ar(k,j)*ti(i,k) return end subroutine to3(cr,ci,rt,it,ar,ai) c C = U* U* U* R implicit none integer*4 i,j,k,l real*8 cr(3,3,3),ci(3,3,3),rt(3,3,3),it(3,3,3), 1ar(3,3),ai(3,3),ti(3,3,3),tr(3,3,3) do 1 i=1,3 do 1 j=1,3 do 1 l=1,3 cr(i,j,l)=0.0d0 ci(i,j,l)=0.0d0 do 1 k=1,3 cr(i,j,l)=cr(i,j,l)+ar(k,i)*rt(k,j,l)+ai(k,i)*it(k,j,l) 1 ci(i,j,l)=ci(i,j,l)-ai(k,i)*rt(k,j,l)+ar(k,i)*it(k,j,l) do 2 i=1,3 do 2 j=1,3 do 2 l=1,3 tr(i,j,l)=0.0d0 ti(i,j,l)=0.0d0 do 2 k=1,3 tr(i,j,l)=tr(i,j,l)+ar(k,j)*cr(i,k,l)+ai(k,j)*ci(i,k,l) 2 ti(i,j,l)=ti(i,j,l)-ai(k,j)*cr(i,k,l)+ar(k,j)*ci(i,k,l) do 3 i=1,3 do 3 j=1,3 do 3 l=1,3 cr(i,j,l)=0.0d0 ci(i,j,l)=0.0d0 do 3 k=1,3 cr(i,j,l)=cr(i,j,l)+ar(k,l)*tr(i,j,k)+ai(k,l)*ti(i,j,k) 3 ci(i,j,l)=ci(i,j,l)-ai(k,l)*tr(i,j,k)+ar(k,l)*ti(i,j,k) return end subroutine toc(mrc,mic,mrt,mit,ar,ai) implicit none integer*4 i,j real*8 mrc(3),mic(3),mrt(3),mit(3),ar(3,3),ai(3,3) do 1 i=1,3 mrc(i)=0.0d0 mic(i)=0.0d0 do 1 j=1,3 c C = Ut* R mrc(i)=mrc(i)+ar(j,i)*mrt(j)+ai(j,i)*mit(j) 1 mic(i)=mic(i)-ai(j,i)*mrt(j)+ar(j,i)*mit(j) return end subroutine tto(r,m,a) implicit none real*8 r(3),m(3),t(3),a(3,3) integer*4 i,j do 1 i=1,3 t(i)=0.0d0 do 1 j=1,3 1 t(i)=t(i)+a(i,j)*r(j) r=t do 3 i=1,3 t(i)=0.0d0 do 3 j=1,3 3 t(i)=t(i)+a(i,j)*m(j) m=t return end subroutine trafo2i(q,a) implicit none integer*4 i,j real*8 q(3,3),t(3,3),a(3,3) c first index, do 11 i=1,3 do 11 j=1,3 11 t(i,j)=a(i,1)*q(1,j)+a(i,2)*q(2,j)+a(i,3)*q(3,j) c second index, do 2 i=1,3 do 2 j=1,3 2 q(i,j)=a(j,1)*t(i,1)+a(j,2)*t(i,2)+a(j,3)*t(i,3) return end subroutine ttQ(r,m,a) implicit none real*8 r(3,3),m(3,3),a(3,3) call trafo2i(r,a) call trafo2i(m,a) return end subroutine trafo3i(q,a) implicit none integer*4 i,jx,kx,j,k real*8 q(3,3,3),t(3,3,3),a(3,3) c first index, do 11 i=1,3 do 11 jx=1,3 do 11 kx=1,3 11 t(i,jx,kx)=a(i,1)*q(1,jx,kx)+a(i,2)*q(2,jx,kx)+a(i,3)*q(3,jx,kx) c second index, do 2 i=1,3 do 2 j=1,3 do 2 kx=1,3 2 q(i, j,kx)=a(j,1)*t(i, 1,kx)+a(j,2)*t(i, 2,kx)+a(j,3)*t(i, 3,kx) c third index, do 3 i=1,3 do 3 j=1,3 do 3 k=1,3 3 t(i, j, k)=a(k,1)*q(i, j, 1)+a(k,2)*q(i, j, 2)+a(k,3)*q(i, j,3) q=t return end subroutine tt3(r,m,a) implicit none real*8 r(3,3,3),m(3,3,3),a(3,3) call trafo3i(r,a) call trafo3i(m,a) return end SUBROUTINE READTTT(ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3) OPEN(15,FILE='FILE.TTT',STATUS='OLD') do 4 i=1,4 4 READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3) CLOSE(15) write(6,*)'ROA tensors read from FILE.TTT' RETURN END c ============================================================ SUBROUTINE TTTT(N,AL,A,G,NAT,ICO,R) c Tensors derivatives C Transforms A and G to local origins (ICO=1) or c common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in in A C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION AL(N,3,3),G(N,3,3),A(N,3,3,3),R(*) real*8,allocatable::X(:,:) allocate(X(3,NAT)) BOHR=0.52917705993d0 DO 1 L=1,NAT DO 1 I=1,3 1 X(I,L)=R(I+3*(L-1))/BOHR C IF(ICO.EQ.1)SI=1.0d0 IF(ICO.EQ.2)SI=-1.0d0 IF(ICO.EQ.3)SI=0.0d0 W=1.0d0 C supposedly the static limit G/w is calculated C DO 2 IA=1,3 DO 2 IB=1,3 DO 2 L=1,NAT DO 2 IE=1,3 II=3*(L-1)+IE call setjk(IB,IG,ID) G(II,IA,IB)=G(II,IA,IB)+ 1SI*0.5d0*W*(X(IG,L)*AL(II,IA,ID)-X(ID,L)*AL(II,IA,IG)) IF(ICO.EQ.3)G(II,IA,IB)=0.0d0 2 CONTINUE IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' C c A IA index is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 II=3*(L-1)+IE A(II,IA,IB,IC)=A(II,IA,IB,IC) 1-1.5d0*SI*(X(IB,L)*AL(II,IA,IC)+X(IC,L)*AL(II,IA,IB)) IF(IB.EQ.IC)A(II,IA,IB,IC)=A(II,IA,IB,IC)+ 1(X(1,L)*AL(II,IA,1)+X(2,L)*AL(II,IA,2)+X(3,L)*AL(II,IA,3))*SI IF(ICO.EQ.3)A(II,IA,IB,IC)=0.0d0 3 CONTINUE IF(ICO.EQ.1)WRITE(6,*)' A transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' A set to zero' RETURN END c c ============================================================ SUBROUTINE TTT2(X,AL,Ga,Aa,Gcom,Ac,natu) implicit none integer*4 natu,A,B,C,G,D,L,IL,I real*8 X(*),AL(natu*3,3,3),Ga(natu*3,3,3),EPS, 1Aa(natu*3,3,3,3),Ac(natu*3,3,3,3),Gcom(natu*3,3,3) c Tensors derivatives C Transforms A and G to common origin c X passed in au C G = static limit G/w DO 1 L=1,natu DO 1 I=1,3 IL=I+3*(L-1) DO 1 A=1,3 DO 1 B=1,3 Gcom(IL,A,B)=Ga(IL,A,B) DO 1 G=1,3 DO 1 D=1,3 1 Gcom(IL,A,B)=Gcom(IL,A,B)-0.5d0*EPS(B,G,D)*X(G+3*(L-1))*AL(IL,A,D) C c A A index is electric DO 2 L=1,natu DO 2 I=1,3 IL=I+3*(L-1) DO 2 A=1,3 DO 2 B=1,3 DO 2 C=1,3 Ac(IL,A,B,C)=Aa(IL,A,B,C) 1+1.5d0*(X(B+3*(L-1))*AL(IL,A,C)+X(C+3*(L-1))*AL(IL,A,B)) 2 IF(B.EQ.C)Ac(IL,A,B,C)=Ac(IL,A,B,C) 1 -X(1+3*(L-1))*AL(IL,A,1) 1 -X(2+3*(L-1))*AL(IL,A,2) 1 -X(3+3*(L-1))*AL(IL,A,3) RETURN END c subroutine tz(alr,ali,gr,gi,ar,ai) implicit none real*8 alr(3,3),ali(3,3),gr(3,3),gi(3,3),ar(3,3,3),ai(3,3,3) alr=0.0d0 ali=0.0d0 gr=0.0d0 gi=0.0d0 ar=0.0d0 ai=0.0d0 return end subroutine setjk(i,j,k) integer*4 i,j,k j=i+1 if(j.gt.3)j=1 k=j+1 if(k.gt.3)k=1 return end subroutine setis(s,is) character*(*) s integer*4 is do 1 is=1,len(s) 1 if(s(is:is).ne.' ')return end subroutine vzd(ur,ui,qr,qi,mr,mi) real*8 ur(*),ui(*),qr(*),qi(*),mr(*),mi(*) integer*4 ix do 13 ix=1,3 ur(ix)=0.0d0 ui(ix)=0.0d0 qr(ix) =0.0d0 qr(ix+3)=0.0d0 qi(ix) =0.0d0 qi(ix+3)=0.0d0 mr(ix)=0.0d0 13 mi(ix)=0.0d0 return end SUBROUTINE writetttd(NAT,ALPHA,G,A) IMPLICIT none INTEGER*4 NAT,I,L,X,J,K,io REAL*8 ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3) io=6 WRITE(io,2000)NAT 2000 FORMAT(' ROA tensors, cartesian derivatives',/,I4,' atoms',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(io,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 X=1,3 1 WRITE(io,2001)L,X,(ALPHA(3*(L-1)+X,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,3g15.7) WRITE(io,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz)') DO 2 I=1,3 WRITE(io,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 X=1,3 2 WRITE(io,2001)L,X,(G(3*(L-1)+X,I,J),J=1,3) WRITE(io,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(io,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 L=1,NAT DO 3 X=1,3 3 WRITE(io,2007)L,X,(A(3*(L-1)+X,I,J,K),K=1,3),L,X,I,J 2007 FORMAT(I5,1H ,I1,3F15.7,' ',4i3) write(io,*) RETURN END subroutine wrx(x,n) implicit none real*8 x(*),BOHR integer*4 n,io,ia,ix io=6 BOHR=0.52917705993d0 do 1 ia=1,n 1 write(io,600)ia,(x(ix+3*(ia-1))*BOHR,ix=1,3) 600 format(i5,3f12.6) return end subroutine wrl(iq,ecm,alr,ali,gr,gi,aar,aai) implicit none integer*4 iq,io,I,J,K real*8 ecm,alr(3,3),ali(3,3),gr(3,3),gi(3,3),aar(3,3,3),aai(3,3,3) io=6 WRITE(io,2003)iq,ecm 2003 FORMAT(i4,f10.2,' cm-1 normal modes derivatives') WRITE(io,*)' Real (alpha) Im (alpha):' DO 11 I=1,3 11 WRITE(io,2001)I,(alr(I,J),J=1,3),(ali(I,J),J=1,3) 2001 FORMAT(I3,6F13.7) WRITE(io,*)' Real (G) Im (G):' DO 21 I=1,3 21 WRITE(io,2001)I,(gr(I,J),J=1,3),(gi(I,J),J=1,3) WRITE(io,*)' Real (A) Im (A):' DO 31 I=1,3 DO 31 J=1,3 31 WRITE(io,2002)I,J,(aar(I,J,K),K=1,3),(aai(I,J,K),K=1,3) 2002 FORMAT(2I3,6F13.7) return end subroutine tcc(n,Cr,Ci,Sr,Si,Eq,Wq) implicit none integer*4 i,j,k,n real*8 Cr(n,n),Ci(n,n),Sr(n,n),Si(n,n),Eq(n),Wq(n),CM real*8,allocatable::tr(:,:),ti(:,:) allocate(tr(n,n),ti(n,n)) CM=219474.630d0 ti=0.0d0 tr=0.0d0 do 1 i=1,n Eq(i)=Wq(i)*CM do 1 j=1,n do 1 k=1,n ti(k,i)=ti(k,i)+Ci(j,i)*Sr(k,j)+Cr(j,i)*Si(k,j) 1 tr(k,i)=tr(k,i)+Cr(j,i)*Sr(k,j)-Ci(j,i)*Si(k,j) Sr=tr Si=ti return end subroutine lookatme(r,i,n,ic) c ic=0: force field c ic=1: Hamiltonian implicit none integer*4 n,ic,im,io,jo,jj,ii real*8 r(n,n),i(n,n),md,od,rii,iii,rij,iij,CM CM=219474.630d0 md=0.0d0 im=0 io=0 jo=0 od=0.0d0 do 1 ii=1,n rii=r(ii,ii) iii=i(ii,ii) if(dabs(rii**2+iii**2).gt.md)then md=rii**2+iii**2 im=ii endif do 1 jj=ii+1,n rij=r(ii,jj) iij=i(ii,jj) if(dabs(rij**2+iij**2).gt.od)then od=rij**2+iij**2 io=ii jo=jj endif 1 continue if(ic.eq.0)then write(6,600)im,im,r(im,im),i(im,im) 600 format('Max FF diag element for i i', 1 2i6,e12.4,' +',e12.4,'i') if(io.eq.0)then write(6,604) else write(6,601)io,jo,r(io,jo),i(io,jo) 601 format('Max FF off-diag element for i j', 1 2i6,e12.4,' +',e12.4,'i') endif else write(6,602)im,im,r(im,im)*CM,i(im,im)*CM 602 format('Max H diag element for i i', 1 2i6,e12.4,' +',e12.4,'i cm-1') if(io.eq.0)then write(6,604) 604 format('Diagonal elements zero') else write(6,603)io,jo,r(io,jo)*CM,i(io,jo)*CM 603 format('Max H off-diag element for i j', 1 2i6,e12.4,' +',e12.4,'i cm-1') endif endif return end subroutine lookate(e,n,ic) c ic=0: w^2 c ic=1: w implicit none integer*4 n,ic,ii real*8 e(n),md,CM CM=219474.630d0 md=0.0d0 do 1 ii=1,n 1 if(e(ii).gt.md)md=e(ii) if(ic.eq.1)md=md*CM write(6,600)md 600 format('Emax: ',f12.4,' cm-1') return end