program aatder implicit none integer*4 nat,I,l,ig98,ndif,NSTEP,N,ia,ix,IC,ltot,ii,ind,jx,k, 1j character*50 s50 character*2 atsym integer*4,allocatable::iz(:),nml(:) real*8,allocatable::r(:,:),AAT(:,:,:),AAT0(:,:,:),S(:,:),E(:), 1stepq(:),rc(:,:),AM(:),aataq(:,:),aptaq(:,:),APT0(:,:,:), 1AATQ(:,:),AAT0Q(:,:),APT0Q(:,:) real*8 ff,ASS,amave,step,ei,ej,f1,bohr,SFACR logical lzmat write(6,6000) 6000 format(/,/, 1' One-purposse program to read AAT derivatives for MVCD',/,/, 2' Input: FILE.OUT',/, 3' S4.OPT ... defined NMDI ndif nml',/, 3' F.INP',/,/, 4' Output: TENAQ.TXT.SCR',/, 5' TENQ.TXT.SCR',/) bohr=0.529177d0 SFACR=0.02342179d0 ndif=0 lzmat=.true. open(10,file='S4.OPT',status='old') 1 read(10,1111,end=99,err=99)s50 IF(s50(1:4).EQ.'ZMAT'.or.s50(2:5).EQ.'ZMAT')THEN READ(10,*)LZMAT ENDIF if(s50(1:4).eq.'NMDI'.or.s50(2:5).eq.'NMDI')then read(10,*)ndif allocate(nml(ndif)) backspace 10 read(10,*)ndif,(nml(i),i=1,ndif) goto 99 endif goto 1 99 close(10) if(ndif.eq.0)call report('ndif = 0') c read S and e, transform into a.u.: open(10,file='F.INP',status='old') read(10,*)nat,nat,nat close(10) N=3*nat allocate(r(nat,3),iz(nat),AAT(nat,3,3),AAT0(nat,3,3),S(N,N), 1E(N),stepq(N),rc(nat,3),AM(N),AATQ(N,3),APT0(nat,3,3), 1aataq(N**2,3),aptaq(2*N**2,3),AAT0Q(N,3),APT0Q(N,3)) call vz(aptaq,2*N**2*3) call vz(aataq,N**2*3) call reads(S,E,N) OPEN(2,FILE='FILE.OUT',status='old') write(6,*)' Looking for geometry' 201 READ(2,1111,END=202)S50 1111 FORMAT(A50) IF((lzmat.and.(s50(19:39).EQ.'Z-Matrix orientation:'.OR. 1 s50(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s50(20:37).EQ.'Input orientation:'.OR. 1 s50(27:44).EQ.'Input orientation:')) 1 .OR. 1 ((.not.lzmat).and. 2 (s50(20:40).EQ.'Standard orientation:'.OR. 2 s50(26:46).EQ.'Standard orientation:')))THEN write(6,1111)S50 ig98=0 if(s50(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s50(27:44).EQ.'Input orientation:'.OR. 1 s50(26:46).EQ.'Standard orientation:')ig98=1 DO 41 I=1,4 41 READ(2,*) do 8 l=1,nat if(ig98.eq.0)then READ(2,*)iz(l),iz(l),(r(l,i),i=1,3) else READ(2,*)iz(l),iz(l),r(l,1),(r(l,i),i=1,3) endif 8 continue goto 203 ENDIF goto 201 202 call report('Geometry not found') 203 continue C c assign masses: do 55 ia=1,nat atsym='??' ASS=amave(iz(IA),atsym) DO 55 I=1,3 55 AM(3*IA-3+I)=ASS NSTEP=1+2*ndif ltot=0 write(6,*)nstep,' harmonic force fields expected' c read-in zero point: c Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0 2 READ(2,1111,END=205)S50 if(s50(2:13).eq.'AAT (total):')then ltot=ltot+1 do 3 ia=1,nat do 3 ix=1,3 3 read(2,*)(AAT0(ia,ix,jx),jx=1,3) goto 206 endif goto 2 205 call report('AAT-0 not found') 206 continue c call vz(stepq,N) c Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0Z0 c start to read the shifted geometries: c DIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIF do 3333 IC=1,2*ndif I=IC if(IC.GT.ndif)I=IC-ndif write(6,*)'Point',IC,'(',I,'):',nml(I),ltot if(IC.le.ndif)then ff=1.0d0 else ff=-1.0d0 endif c c check, that we are at the right point, eventually extract step: step=0.0d0 call chkgeos(r,nat,lzmat,step,N,S,rc,nml(I),AM,E) stepq(nml(I))=step c c read-in current values: 6 READ(2,1111,END=208)S50 if(s50(2:13).eq.'AAT (total):')then ltot=ltot+1 do 7 ia=1,nat do 7 ix=1,3 7 read(2,*)(AAT(ia,ix,jx),jx=1,3) goto 207 endif goto 6 208 call report('AAT not found') 207 continue c transform dipole derivatives into normal modes: call trafopq(AAT,AATQ,S,N) DO 543 II=1,N DO 543 K=1,3 ind=nml(I)+N*(II-1) 543 aataq(ind,K)=aataq(ind,K)+AATQ(II,K)*ff c aataq - nml(I) index goes to coordinate, II is momentum c 3333 continue c DIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIF close(2) C if(stepq(nml(I)).lt.1.0d-9)call report('Diff. step too small!') c DO 18 I=ndif,1,-1 ei=sqrt(E(nml(I))) f1=0.5d0* SFACR*bohr/stepq(nml(I)) DO 18 J=ndif,1,-1 ej=sqrt(E(nml(J))) ind=nml(I)+N*(nml(J)-1) do 18 k=1,3 18 aataq(ind,k)=aataq(ind,k)/ei*ej*f1 c aataq = i-index coordinate, J -momentum call wrtenq(ndif,N,aptaq,aataq,nml) open(10,file='FILE.TEN',status='old') read(10,*) do 9 ia=1,nat do 9 ix=1,3 9 read(10,*)(APT0(ia,ix,jx),jx=1,3) do 10 ia=1,nat do 10 ix=1,3 10 read(10,*)(AAT0(ia,ix,jx),jx=1,3) close(10) call trafopq(AAT0,AAT0Q,S,N) call trafopq(APT0,APT0Q,S,N) DO 161 I=ndif,1,-1 ei=sqrt(E(nml(I))) do 161 j=1,3 APT0Q(nml(I),j)=APT0Q(nml(I),j)/ei 161 AAT0Q(nml(I),j)=AAT0Q(nml(I),j)*ei call wrten1q(ndif,APT0Q,AAT0Q,nml,N) end c SUBROUTINE READS(S,E,N) implicit none integer*4 N,NM,nat,I,IA,K,iq real*8 S(N,N),E(*),cm,SFACR cm=219470.0d0 SFACR=0.02342179d0 OPEN(7,FILE='F.INP') READ(7,*)NM,N,nat DO 2 I=1,N E(I)=0.0d0 DO 2 IA=1,N 2 S(IA,I)=0.0d0 DO 1 I=1,nat 1 READ(7,*) READ(7,*) DO 3 IA=1,NAT DO 3 I=1,NM 3 READ(7,*)(S(3*(IA-1)+K,N-I+1),K=1,2),(S(3*(IA-1)+K,N-I+1),K=1,3) READ(7,*) READ(7,*)(E(N-I+1),I=1,NM) CLOSE(7) do 9 iq=1,nm E(N-iq+1)=E(N-iq+1)/cm do 9 i=1,nat do 9 k=1,3 9 S(3*(i-1)+k,N-iq+1)=S(3*(i-1)+k,N-iq+1)*SFACR write(6,*)'S-matrix has been read from F.INP, changed to au' RETURN END subroutine trafopq(AAT,AATQ,S,N) implicit none integer*4 N,i,j,a,b real*8 AAT(N/3,3,3),AATQ(N,3),S(N,N) DO 1 i=1,3 DO 1 a=1,N AATQ(a,i)=0.0d0 DO 1 b=1,N/3 DO 1 j=1,3 1 AATQ(a,i)=AATQ(a,i)+S(3*(b-1)+j,a)*AAT(b,j,i) return end subroutine chkgeos(r,nat,lzmat,step,N,S,rc,ireq,AM,E) implicit none integer*4 nat,ig98,I,l,ic,ix,iq,ii,ireq,N real*8 r(nat,3),step,rc(nat,3),dq,S(N,N),CM,E(*), 1SFACR,bohr,geotol,AM(*) character*50 S50 logical lzmat c CM=219470.0d0 SFACR=0.02342179d0 bohr=0.529177d0 geotol=1.0d-4 c 201 READ(2,1111,END=202)S50 1111 FORMAT(A50) IF((lzmat.and.(s50(19:39).EQ.'Z-Matrix orientation:'.OR. 1 s50(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s50(20:37).EQ.'Input orientation:'.OR. 1 s50(27:44).EQ.'Input orientation:')) 1 .OR. 1 ((.not.lzmat).and. 2 (s50(20:40).EQ.'Standard orientation:'.OR. 2 s50(26:46).EQ.'Standard orientation:')))THEN ig98=0 if(s50(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s50(27:44).EQ.'Input orientation:'.OR. 1 s50(26:46).EQ.'Standard orientation:')ig98=1 DO 41 I=1,4 41 READ(2,*) ic=0 c c read current geometry DO 1 l=1,nat if(ig98.eq.0)then READ(2,*)I,I,(rc(l,ix),ix=1,3) else if(ig98.eq.-1)then READ(2,*)I,(rc(l,ix),ix=1,3) else READ(2,*)I,I,I,(rc(l,ix),ix=1,3) endif endif 1 continue c c Cartesian displacements from equilibrium DO 2 l=1,nat DO 2 ix=1,3 c write(6,*)l,ix,rc(l,ix),r(l,ix) 2 rc(l,ix)=rc(l,ix)-r(l,ix) c c normal mode displacements do 3 iq=1,3*nat dq=0.0d0 do 4 l=1,nat do 4 ix=1,3 ii=3*(l-1)+ix c write(6,*)ii,iq,ix,l,AM(ii),S(ii,iq),rc(l,ix) 4 dq=dq+AM(ii)*S(ii,iq)*rc(l,ix) c c adapt to the strange pmz units dq=dq/SFACR: dq=dq/SFACR c stop if(dabs(dq).gt.geotol)then write(6,*)nat write(6,6000)iq,e(iq)*CM,dq,dq/(SFACR*bohr) 6000 format('mode',i4,f10.2,' cm-1',' dQ =',f10.4,' au:',f10.4) if(ireq.ne.iq)call report('errorneous mode numbering') ic=ic+1 call chksti(dq,step) endif 3 continue c if(ic.ne.1)then write(3,*)ic call report('deffective differentiation') endif return ENDIF goto 201 202 call report('Geometry not found') end subroutine chksti(d,step) real*8 d,step if(step.lt.1.0d-9)then step=dabs(d) write(3,3000)step write(6,3000)step 3000 format('Step set on the fly as ',F10.5) endif return end c ============================================================ function amave(iz,sym) c atomic mass of element with atomiz number iz c approximate - to be compatible with new2 implicit none integer*4 iz,MENDELEV PARAMETER (MENDELEV=89) CHARACTER*2 atsy(MENDELEV),sym data atsy/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne', 3'Na','Mg','Al','Si',' P',' S','Cl','Ar', 4' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn', 4 'Ga','Ge','As','Se','Br','Kr', 5'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd', 5 'In','Sn','Sb','Te',' I','Xe', 6'Cs','Ba','La', 6 'Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho', 6 'Er','Tm','Yb','Lu', 6'Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg', 6 'Tl','Pb','Bi','Po','At','Rn', 7'Fr','Ra','Ac'/ real*8 amave,amas(MENDELEV) data amas/1.007825d0,4.003, 2 6.941, 9.012, 10.810,12.0d0,14.003074d0,15.999d0, 2 18.9984d0, 20.179, 3 22.990,24.305, 26.981,28.086,30.9737d0, 3 31.972d0,34.968852d0,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ c if(iz.lt.0.or.iz.gt.MENDELEV)then write(6,*)iz call report('Unknown atomic number') else amave=amas(iz) sym=atsy(iz) return endif end subroutine report(s) character*(*) s write(6,*)s stop end subroutine vz(v,N) real*8 v(*) integer*4 i,n do 1 i=1,n 1 v(i)=0.0d0 return end subroutine wrtenq(ndif,N,aptaq,aataq,nml) implicit none integer*4 I,ndif,nml(*),J,ind1,ind2,k,N real*8 aptaq(2*N**2,3),aataq(N**2,3) open(9,file='TENAQ.TXT.SCR') write(9,*)ndif write(9,9002) 9002 format(' APT second derivatives',/, 1' Iq Jq ix PIqJq PJqIq PIqIqJq') DO 18 I=ndif,1,-1 DO 18 J=I,1,-1 ind1=nml(I)+N*(nml(J)-1) ind2=nml(J)+N*(nml(I)-1) do 18 k=1,3 18 write(9,9003)nml(I),nml(J),k,aptaq(ind1,k),aptaq(ind2,k), 1aptaq(N**2+ind1,k) 9003 format(3i5,3g15.5) write(9,9004) 9004 format(' AAT second derivatives',/, 1' Iq Jp ix PIqJp PJpIq') DO 1 I=ndif,1,-1 DO 1 J=I,1,-1 ind1=nml(I)+N*(nml(J)-1) ind2=nml(J)+N*(nml(I)-1) do 1 k=1,3 1 write(9,9003)nml(I),nml(J),k,aataq(ind1,k),aataq(ind2,k) close(9) write(6,*)'Second dipole derivatives written to TENAQ.TXT.SCR' return end subroutine wrten1q(ndif,APT0Q,AAT0Q,nml,N) implicit none integer*4 I,j,ndif,nml(*),N real*8 APT0Q(N,3),AAT0Q(N,3) open(9,file='TENQ.TXT.SCR') write(9,9001)0.0d0,0.0d0,0.0d0 9001 format(3f15.7,' dipole (au)',/, 1' Normal mode derivatives of dipole moment num num-second anal') write(9,*)ndif DO 16 I=ndif,1,-1 do 16 j=1,3 16 write(9,9000)I,nml(I),j,APT0Q(nml(I),j),0.0d0,APT0Q(nml(I),j) 9000 format(3i5,3f12.6) write(9,9022) 9022 format(' AAT:') DO 17 I=ndif,1,-1 do 17 j=1,3 17 write(9,9000)I,nml(I),j,AAT0Q(nml(I),j) close(9) write(6,*)'APT AAT in TENQ.TXT.SCR' return end