PROGRAM analFINP IMPLICIT REAL*8(A-H,O-Z) integer,allocatable::q(:),bt(:,:),nt(:),a(:,:) real*8,allocatable::r(:,:),s(:,:),u(:,:),e(:) open(10,file='FILE.X',status='old') read(10,*) read(10,*)nat allocate(q(Nat),r(Nat,3),bt(Nat,7),nt(Nat),a(Nat,8)) do 101 i=1,NAT read(10,*)q(i),(r(i,j),j=1,3),(bt(i,j),j=1,7) nt(i)=0 do 101 j=1,7 101 if(bt(i,j).ne.0)nt(i)=nt(i)+1 close(10) write(6,*)'FILE.X read,',nat,' atoms' c identify amides do 104 na=1,nat do 104 ii=1,8 104 a(na,ii)=0 na=0 do 102 ia=1,nat if(q(ia).eq.8.and.nt(ia).eq.1)then i1=bt(ia,1) if(q(i1).eq.6.and.nt(i1).eq.3)then il=0 in=0 do 103 i2=1,nt(i1) ii=bt(i1,i2) if(q(ii).eq.7)in=ii 103 if(q(ii).eq.6)il=ii if(in.ne.0)then c a1 a2 a3 a4 a5 a6 a7 a8 c C -- N -- aC -- C ( = O) -- N -- aC -- C c i11 in1 il i1 ia in il2 i22 na=na+1 in1=0 do 105 ii=1,nt(il) jj=bt(il,ii) 105 if(q(jj).eq.7.and.nt(jj).eq.3)in1=jj i11=0 if(in1.ne.0)then do 106 ii=1,nt(in1) jj=bt(in1,ii) if(q(jj).eq.6.and.jj.ne.il)then do 110 kk=1,nt(jj) k=bt(jj,kk) 110 if(q(k).eq.8.and.nt(k).eq.1)i11=jj endif 106 continue endif il2=0 if(in.ne.0)then do 107 ii=1,nt(in) jj=bt(in,ii) if(q(jj).eq.6.and.jj.ne.i1)then do 109 kk=1,nt(jj) k1=bt(jj,kk) 109 if(q(k1).eq.6.and.nt(k1).eq.3)il2=jj endif 107 continue endif i22=0 if(il2.ne.0)then do 108 ii=1,nt(il2) jj=bt(il2,ii) 108 if(q(jj).eq.6.and.nt(jj).eq.3)i22=jj endif a(na,1)=i11 a(na,2)=in1 a(na,3)=il a(na,4)=i1 a(na,5)=ia a(na,6)=in a(na,7)=il2 a(na,8)=i22 write(6,800)na,(a(na,j),j=1,8) 800 format(i4,8i5) endif endif endif 102 continue write(6,*)na,' amides' open(4,file='F.INP',status='old') read(4,*)M,N,Nat allocate(s(N,M),u(M,6),e(M)) write(6,*)' F.INP: ',M,' modes, ',Nat,' atoms' do 1 i=1,NAT 1 read(4,*)q(i),(r(i,j),j=1,3) read(4,*) DO 2 I=1,NAT DO 2 J=1,M 2 read(4,*)(s(3*(i-1)+ix,j),ix=1,2),(s(3*(i-1)+ix,j),ix=1,3) read(4,*) READ(4,4000)(E(I),I=1,M) 4000 FORMAT(6F11.6) do 3 i=1,M 3 read(4,*)(u(i,j),j=1,6) close(4) write(6,*)'F.INP read',nat,' atoms' do 4 i=1,M c amide with largest shift i1=1 d1=0.0 do 5 ia=1,na ic=a(ia,4) io=a(ia,5) sx=s(3*(ic-1)+1,i)-s(3*(io-1)+1,i) sy=s(3*(ic-1)+2,i)-s(3*(io-1)+2,i) sz=s(3*(ic-1)+3,i)-s(3*(io-1)+3,i) rx=r(ic,1)-r(io,1) ry=r(ic,2)-r(io,2) rz=r(ic,3)-r(io,3) da2=sx*rx+sy*ry+sz*rz d2=rx*rx+ry*ry+rz*rz da=da2/sqrt(d2) if(da.gt.d1)then d1=da i1=ia endif 5 continue if(d1.gt.0.061.and.e(i).gt.1650.)then i8=0 do 111 ii=1,8 111 if(a(i1,ii).ne.0)i8=i8+1 if(i8.eq.8)then w=tor(r,nat,a(i1,3),a(i1,4),a(i1,6),a(i1,7)) p=tor(r,nat,a(i1,2),a(i1,3),a(i1,4),a(i1,6)) f=tor(r,nat,a(i1,1),a(i1,2),a(i1,3),a(i1,4)) g=tor(r,nat,a(i1,4),a(i1,6),a(i1,7),a(i1,8)) write(6,605)i,e(i),i1,d1,w,p,f,g 605 format(' Amide mode',i6,f10.2,' cm-1, group ',i6,f10.4,4f8.1) endif endif 4 continue end function di(r,Nat,i,j) dimension r(Nat,3) d2=(r(i,1)-r(j,1))**2+(r(i,1)-r(j,1))**2+(r(i,1)-r(j,1))**2 di=sqrt(d2) return end function tor(x,nat,i1,i2,i3,i4) IMPLICIT REAL*8(A-H,O-Z) dimension x(nat,3) DIMENSION XX(3),YY(3),A(3),P(3),E0(3),B0(3), 1b(3),c(3),v1(3),v2(3),AP(3) logical lex,ldistr,lper PI=3.141592653589 C do 23 ix=1,3 A(ix) =X(i3,ix)-X(i4,ix) P(ix) =X(i1,ix)-X(i2,ix) 23 E0(ix)=X(i2,ix)-X(i3,ix) call vp(B0,E0,A) SSB0B0=sp(B0,B0) SSE0E0=sp(E0,E0) DO 8 KK=1,3 B0(KK)=B0(KK)/DSQRT(SSB0B0) 8 E0(KK)=E0(KK)/DSQRT(SSE0E0) SSPP=sp(P,P) SSAA=sp(A,A) SSAP=sp(A,P) AD=DSQRT(SSPP) SSAE0=sp(A,E0) SSPE0=sp(P,E0) DO 9 KK=1,3 XX(KK)=A(KK)-E0(KK)*SSAE0 9 YY(KK)=P(KK)-E0(KK)*SSPE0 SSXXXX=sp(XX,XX) SSYYYY=sp(YY,YY) SSXXYY=sp(XX,YY) call vp(AP,XX,YY) SSE0AP=sp(AP,E0) SIGN=1.0D0 IF (SSE0AP.NE.0.0D0)SIGN=-SSE0AP/ABS(SSE0AP) CB=0.0D0 IF (AD.NE.0.0D0)CB=-SSPE0/AD CA=0.0D0 IF (SSXXXX.NE.0.0D0.AND.SSYYYY.NE.0.0D0) 1CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY) ANGLE=0.0D0 IF (CA.LE.-1.0D0) ANGLE=PI IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA) ANGLE=180.0D0* ANGLE/PI ANGLE=SIGN*ANGLE if(ANGLE.lt.-180.0d0)ANGLE= 360.0d0+ANGLE if(ANGLE.gt.+180.0d0)ANGLE=-360.0d0+ANGLE tor=ANGLE return end subroutine report(s) character*(*) s write(6,*)s stop end function sp(a,b) IMPLICIT none real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine vp(A,X,Y) implicit none real*8 A(*),X(*),Y(*) A(1)=X(2)*Y(3)-X(3)*Y(2) A(2)=X(3)*Y(1)-X(1)*Y(3) A(3)=X(1)*Y(2)-X(2)*Y(1) return end