program ped implicit none integer*4 NOB,NA,nat,I,nn,ia,ix,a,ii,npe,imax,j,NQ,iargc real*8 amax,anorm,amf real*8,allocatable::B(:,:),X(:),S(:,:),FF(:,:),SB(:,:),w(:), 1AA(:,:),AMODE(:),PE(:),AAT(:,:) integer*4,allocatable::MODE(:),IPE(:),iz(:) character*1 fok write(6,6000) 6000 format(/,' Potential energy distribution / mixed calculation',/, 1 /,' Input: B.MAT Cartesian-Internal ', 1 /,' transformation matrix', 1 /,' INTY.FC internal force field', 1 /,' F.INP S-matrix',/, 1 /,' Output: POT_all.LST PED - all', 1 /,' POT.EN PED - selected',/,/) fok='Y' if(iargc().eq.1)call getarg(1,fok) if(fok.eq.'y')fok='Y' OPEN(9,FILE='B.MAT',status='old') read(9,*)NOB read(9,*)NA allocate(B(NOB,NA)) B=0.0d0 read(9,*) do 1 I=1,NOB read(9,*) read(9,*)nn do 1 ii=1,nn 1 read(9,*)ia,ix,B(I,ix+3*(ia-1)) CLOSE(9) nat=NA/3 write(6,6001)nat,NOB 6001 format(' B.MAT read,',i6,' atoms,',i6,' internal coordinates') allocate(X((NOB*(NOB+1))/2),FF(NOB,NOB)) call riff(NOB,X,FF,fok) allocate(S(NA,NA),w(NA),iz(nat)) call reads(NQ,nat,S,w,iz) allocate(SB(NQ,NOB)) do 2 i=1,NQ do 2 j=1,NOB SB(i,j)=0.0d0 do 2 a=1,NA ia=(a+2)/3 2 SB(i,J)=SB(i,J)+S(a,i)*B(J,a)*amf(iz(ia)) allocate(AA(NOB,NQ),AAT(NOB,NQ)) do 3 i=1,NQ do 3 J=1,NOB if(dabs(w(i)).gt.1.0d0)then AA(J,i)=SB(i,J)**2*dabs(FF(J,J))/w(i) else AA(J,i)=SB(i,J)**2*dabs(FF(J,J)) endif 3 continue do 6 i=1,NQ ANORM=0.0d0 do 7 J=1,NOB 7 ANORM=ANORM+AA(J,i) if(dabs(ANORM).gt.1.0d-14)then do 61 J=1,NOB 61 AA(J,i)=AA(J,i)/ANORM endif 6 continue OPEN(17,FILE='POT_all.LST') do 4 i=1,NQ write(17,1717)i,w(i) 1717 format(i6,f10.2,$) do 5 J=1,NOB 5 write(17,17171)AA(J,i) 17171 format(f10.6,$) 4 write(17,*) close(17) write(6,*)'POT_all.LST written' allocate(AMODE(NOB),MODE(NOB),PE(NOB),IPE(NOB)) OPEN(19,FILE='POT.EN') WRITE(19,*)' POTENTIAL ENERGY DISTRIBUTION [in %]' WRITE(19,*)' MODE coordinate(PED) (for PED > 5%)' WRITE(19,1920) 1920 FORMAT(80(1H-)) DO 454 J=1,NOB AMODE(J)=0.0D0 454 MODE(J)=0 c save for the last analysis: AAT=AA c c Loop over normal modes: DO 451 J=1,NQ NPE=0 c c Loop over coordinates DO 453 II=1,NOB AMAX=AA(1,J) IMAX=1 c c Second loop over coordinates - find maximum and zero after DO 452 I=2,NOB IF(AA(I,J).GT.AMAX)THEN AMAX=AA(I,J) IMAX=I ENDIF 452 CONTINUE AA(IMAX,J)=0.0D0 AMAX=AMAX*100.0D0 c IF (AMAX.GE.5.0D0)THEN NPE=NPE+1 PE(NPE)=AMAX IPE(NPE)=IMAX ENDIF c c record the absolute maximum:(mode J, coord IMAX) IF(AMODE(IMAX).LT.AMAX)THEN MODE(IMAX)=J AMODE(IMAX)=AMAX ENDIF 453 CONTINUE NPE=MIN(NPE,5) if(PE(1).gt.99.9d0)PE(1)=99.9d0 IF(NPE.EQ.1)WRITE(19,1922)J,w(J),IPE(1),PE(1) IF(NPE.EQ.2)WRITE(19,1922)J,w(J),IPE(1),PE(1), 1IPE(2),PE(2) IF(NPE.EQ.3)WRITE(19,1922)J,w(J),IPE(1),PE(1), 1IPE(2),PE(2),IPE(3),PE(3) IF(NPE.EQ.4)WRITE(19,1922)J,w(J),IPE(1),PE(1), 1IPE(2),PE(2),IPE(3),PE(3),IPE(4),PE(4) IF(NPE.EQ.5)WRITE(19,1922)J,w(J),IPE(1),PE(1), 1IPE(2),PE(2),IPE(3),PE(3),IPE(4),PE(4),IPE(5),PE(5) 1922 FORMAT(I5,F9.1, 5(I5,'(',F4.1,')') ) 451 CONTINUE WRITE(19,1920) DO 455 I=1,NOB MODE(I)=1 do 8 J=1,NQ c take the saved matrix: 8 if(dabs(AAT(I,J)).gt.dabs(AAT(I,MODE(I))))MODE(I)=J 455 WRITE(19,1919)I,MODE(I),w(MODE(I)),AAT(I,MODE(I))*100.0d0 1919 FORMAT(' Coordinate',I4,' is most present in mode',I4, 1f10.2, 'cm-1 (by',F5.1,'%).') WRITE(19,1920) CLOSE(19) write(6,*)'POT.EN written' end subroutine reads(NQ,NOAT,S,w,iz) implicit none integer*4 NQ,NOAT,K,I,J,iz(*) real*8 S(3*NOAT,3*NOAT),w(*) OPEN(17,FILE='F.INP',STATUS='OLD') READ(17,*)NQ DO 3 K=1,NOAT 3 READ(17,*)iz(K) READ(17,*) DO 4 K=1,NOAT DO 4 I=1,NQ 4 READ (17,*)(S(J+3*(K-1),I),J=1,2),(S(J+3*(K-1),I),J=1,3) READ(17,*) READ(17,*)(w(I),I=1,NQ) close(17) j=0 do 1 i=1,NQ 1 if(w(i).gt.1.0d0)j=j+1 WRITE(*,*)' S-matrix read in,',NQ,' modes,',j,' are positive' return end subroutine riff(NOB,X,FF,fok) implicit none integer*4 N,N1,N3,LN,J,I,K,NOB real*8 X(*),FF(NOB,NOB) character*1 fok if(fok.eq.'Y') then OPEN(8,FILE='INTY.FC') read(8,*) N=0 N1=1 2 N3=N1+4 IF(N3.GT.NOB)N3=NOB read(8,*) DO 3 LN=N1,NOB read(8,*)X((LN*(LN-1))/2+N1),(X((LN*(LN-1))/2+J),J=N1,MIN(LN,N3)) 3 N=N+MIN(LN,N3)-N1+1 N1=N1+5 IF(N3.LT.NOB)GOTO 2 close(8) K=0 DO 4 I=1,NOB DO 4 J=1,I K=K+1 FF(I,J)=X(K) 4 FF(J,I)=X(K) write(6,*)'INTY.FC read' else FF=0.0d0 DO 1 I=1,NOB 1 FF(I,I)=1.0d0 write(6,*)'INTY.FC set to unity for distance distribution.' endif return end function amf(z) implicit none real*8 amf integer*4 z,MENDELEV PARAMETER (MENDELEV=89) real*4 AM(MENDELEV) data AM/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,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/ if(z.lt.1.or.z.gt.MENDELEV)then write(6,*)z,'invalid atomic number' stop endif amf=dble(AM(z)) return end