program nmaai3 implicit none integer*4 NAT,ia,io,ix,nao,NATS,IERR,N3,NM,N,i,j,is,jx,i1,i2,no, 1m1,m2,mm,im,i1i,i1j,i2i,i2j,mi,oi,oj,j1,j2,mj,nn integer*4,allocatable::ND(:),IBONDSB(:,:),NBOB(:),INDBIG(:,:), 1NDS(:),IBONDS(:,:),NBO(:) real*8,allocatable::R(:,:),XS(:,:),XB(:,:),RS(:,:),S(:,:),E(:), 1F(:,:),SB(:,:),EB(:),fi(:,:),t(:,:),AAT(:,:,:),P(:,:,:),AAC(:,:,:) real*8 A(3,3),sx,sy,sz,CM,wi,wj,hij,u(3),m(3),DS,RD,c(3) character*8 s8,si,sj,sii,sjj logical lex write(6,6009) 6009 format(/,' Gets interaction for AI modes based on NMA ',/, 1' in a tripeptide',/, 2' Input: OVER.TXT',/, 3' BIG.X',/, 5' BIG.FC',/, 4' BIG.TEN (optional)',/, 4' SMALL.X',/, 4' SMALLF.INP',/) m1=8 m2=8 no=2 CM=219474.0d0 if(iargc().eq.2)then call getarg(1,s8) read(s8,*)m1 read(s8,*)m2 call getarg(2,s8) read(s8,*)no endif if(iargc().eq.3)then call getarg(1,s8) read(s8,*)m1 call getarg(2,s8) read(s8,*)m2 call getarg(3,s8) read(s8,*)no endif NAT=1 allocate(R(3,1),ND(1),IBONDSB(1,1),NBOB(1)) CALL READRX('BIG.X',NAT,R,ND,IBONDSB,NBOB,0) deallocate(R,ND,IBONDSB,NBOB) allocate(R(3,NAT),ND(NAT),IBONDSB(4,NAT),NBOB(NAT)) CALL READRX('BIG.X',NAT,R,ND,IBONDSB,NBOB,1) WRITE(6,*)NAT, ' atoms in the big molecule' NATS=1 allocate(RS(3,1),NDS(1),IBONDS(1,1),NBO(1)) CALL READRX('SMALL.X',NATS,RS,NDS,IBONDS,NBO,0) deallocate(RS,NDS,IBONDS,NBO) allocate(RS(3,NATS),NDS(NATS),IBONDS(4,NATS),NBO(NATS)) CALL READRX('SMALL.X',NATS,RS,NDS,IBONDS,NBO,1) WRITE(6,*)NATS, ' atoms in the small molecule' allocate(INDBIG(no,NAT)) call RDOVER(INDBIG,NAT,NDS,ND,no) N3=3*NATS allocate(S(N3,N3),E(N3)) call READS('SMALLF.INP',N3,S,E,NM) write(6,888)m1,E(m1),m2,E(m2),no 888 format(' Modes ',i4,'(',f12.2,') - ',i4,'(',f12.2,'(;', 1i3,' overlaps') N=3*NAT allocate(F(N,N)) call READFF('BIG.FC',N,F) mm=(m2-m1+1)*no allocate(SB(N,mm),EB(mm)) do 5 i=1,mm do 5 j=1,N 5 SB(j,i)=0.0d0 inquire(file='BIG.TEN',exist=lex) if(lex)then allocate(P(NAT,3,3),AAT(NAT,3,3),AAC(NAT,3,3)) CALL READTEN('BIG.TEN',NAT,P,AAT) endif do 1 io=1,no write(6,*)' Overlap ',io nao=0 do 2 ia=1,NAT 2 if(INDBIG(io,ia).ne.0)nao=nao+1 write(6,*)nao,' overlaped atoms' allocate(XS(nao,3),XB(nao,3)) nao=0 do 3 ia=1,NAT if(INDBIG(io,ia).ne.0)then nao=nao+1 do 31 ix=1,3 XS(nao,ix)=RS(ix,INDBIG(io,ia)) 31 XB(nao,ix)=R( ix,ia) endif 3 continue call center(XS,nao) call center(XB,nao) A(1,1)=0.0d0 call DOMATRIX(IERR,A,XS,XB,nao) deallocate(XS,XB) c c NMA S-matrix rotated to amide group in the dippetide, c and with the dipeptide atoms, amide I = # 8 in NMA do 4 ia=1,NAT is=INDBIG(io,ia) if(is.ne.0)then do 11 im=m1,m2 j=io+no*(im-m1) sx=A(1,1)*S(3*is-2,im)+A(1,2)*S(3*is-1,im)+A(1,3)*S(3*is,im) sy=A(2,1)*S(3*is-2,im)+A(2,2)*S(3*is-1,im)+A(2,3)*S(3*is,im) sz=A(3,1)*S(3*is-2,im)+A(3,2)*S(3*is-1,im)+A(3,3)*S(3*is,im) SB(3*ia-2,j)=sx SB(3*ia-1,j)=sy 11 SB(3*ia ,j)=sz endif 4 continue do 12 im=m1,m2 12 EB(io+no*(im-m1))=E(im) if(lex)then write(si,608)io call i1i2(i1,i2,si) do 15 im=m1,m2 j=io+no*(im-m1) write(sii,608)im call i1i2(i1i,i2i,sii) c electric and magnetic dipole moments: call dd(u,m,P,AAT,NAT,nm,E,im,DS,RD,SB,io,j) c amide center: do 14 ix=1,3 c(ix)=0.0d0 nn=0 do 141 ia=1,NAT is=INDBIG(io,ia) if(is.ne.0)then nn=nn+1 c(ix)=c(ix)+R(ix,ia) endif 141 continue 14 c(ix)=c(ix)/dble(nn) c put AAT to amide center: call PO(NAT,P,AAT,AAC,c) write(6,*)'In amide center:' call dd(u,m,P,AAC,NAT,nm,E,im,DS,RD,SB,io,j) open(7,file='r'//si(i1:i2)//'.'//sii(i1i:i2i)) write(7,702)RD 702 format(e11.4) 15 close(7) endif write(6,*)' A:' 1 write(6,618)A 618 format(3f12.6) call wrf('CONTROL.F.INP',SB,R,NAT,mm,EB,ND) allocate(t(mm,N)) do 61 i =1,mm do 61 jx=1,N t(i,jx)=0.0d0 do 61 ix=1,N 61 t(i,jx)=t(i,jx)+SB(ix,i)*F(ix,jx) allocate(fi(mm,mm)) do 62 i =1,mm do 62 j =1,mm fi(i,j)=0.0d0 do 62 jx=1,N 62 fi(i,j)=fi(i,j)+SB(jx,j)*t(i,jx) open(71,file='VIJ.TXT') do 63 i=1,mm wi=dsqrt(fi(i,i)) mi=m1+(i+no-1)/no-1 write(sii,608)mi call i1i2(i1,i2,sii) oi=i-(mi-m1)*no write(si,608)oi 608 format(i8) call i1i2(i1i,i2i,si) do 63 j=i,mm wj=dsqrt(fi(j,j)) mj=m1+(j+no-1)/no-1 write(sjj,608)mj call i1i2(j1,j2,sjj) oj=j-(mj-m1)*no write(sj,608)oj call i1i2(i1j,i2j,sj) if(i.eq.j)then hij=wi*CM else hij=fi(i,j)/dsqrt(wi*wj)*CM/2.0d0 endif write(6,607)i,j,hij 607 format(' H',i2,i2,': ',F12.2) if(m1.eq.m2)then c one mode: open(7,file='v'//si(i1i:i2i)//sj(i1j:i2j)//'.'//sii(i1:i2)) else c more modes: open(7,file='v'//si(i1i:i2i)//sj(i1j:i2j) 1 //'.'//sii(i1:i2)//'.'//sjj(j1:j2)) endif write(7,700)hij 700 format(f14.4) write(71,701)i,j,hij 701 format(2i4,f14.4) 63 close(7) end subroutine i1i2(i1,i2,s) implicit none character(*) s integer*4 i1,i2 do 1 i1=1,len(s) 1 if(s(i1:i1).ne.' ')goto 2 2 do 3 i2=len(s),1,-1 3 if(s(i2:i2).ne.' ')goto 4 4 return end subroutine wrf(f,S,R,NAT,NM,EB,ND) implicit none character(*) f integer*4 NAT,NM,ND(*),I,ix,IU,J,k real*8 S(3*NAT,NM),EB(*),R(3,NAT),SFACR SFACR=0.02342179d0 OPEN(34,FILE=f) WRITE(34,10)NM,3*NAT,NAT 10 FORMAT(3I7) DO 3 I=1,NAT 3 WRITE(34,11)ND(I),(R(ix,I),ix=1,3) 11 FORMAT(I7,3F12.6) WRITE(34,14) 14 FORMAT(' Atom Mode X-disp. Y-disp. Z-disp.') DO 1 I=1,NAT DO 1 J=1,NM IU=3*(I-1) 1 WRITE(34,15)I,J,(S(IU+k,J)/SFACR,k=1,3) 15 FORMAT(2I7,3F11.6) WRITE(34,11)NM WRITE(34,13)(EB(I),I=1,NM) 13 format(6F11.3) CLOSE(34) RETURN END SUBROUTINE READRX(f,NAT,R,ND,IBONDS,NB,ic) IMPLICIT none integer*4 NAT,ND(*),IBONDS(4,NAT),NB(*),IC,i,IB REAL*8 R(3,NAT) character*(*) f OPEN(2,FILE=f,STATUS='OLD') READ(2,*) read(2,*)NAT if(ic.eq.0)goto 99 do 1 i=1,NAT READ(2,*)ND(i),R(1,i),R(2,i),R(3,i),(IBONDS(IB,i),IB=1,4) NB(i)=0 DO 1 IB=1,4 1 IF(IBONDS(IB,i).NE.0)NB(i)=IB 99 CLOSE(2) RETURN END SUBROUTINE RDOVER(INDBIG,NAT,NDS,ND,no) implicit none integer*4 NAT,no,INDBIG(no,NAT),I,IAB,NDS(*),ND(*) open(9,file='OVER.TXT',status='old') do 1 I=1,no READ(9,90)(INDBIG(I,IAB),IAB=1,NAT) 90 format(20i3) READ(9,*)(INDBIG(I,IAB),IAB=1,NAT) DO 1 IAB=1,NAT IF(INDBIG(I,IAB).NE.0)then IF(ND(IAB).NE.NDS(ABS(INDBIG(I,IAB))))THEN WRITE(6,3000)I,IAB,ND(IAB),NDS(ABS(INDBIG(I,IAB))) 3000 FORMAT('Overlap ',I2,', atom ',I3,': ',2I4) CALL REPORT(' Different atomic types !') ENDIF ENDIF 1 CONTINUE close(9) write(6,*)'overlap read' return end SUBROUTINE REPORT(RS) CHARACTER*(*) RS WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(6,*)RS WRITE(6,3001) 3001 FORMAT(80(1H*),/,/,'PROGRAM STOPPED') STOP END SUBROUTINE DOMATRIX(IERR,A,XS,XB,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 NAT,ITER real*8 ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3), 1A(3,3),XS(NAT,3),XB(NAT,3),RTOL IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG,A,XS,XB,NAT) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG,A,XS,XB,NAT) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG,A,XS,XB,NAT) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG,A,XS,XB,NAT) FTOL=0.0000001d0 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-11)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)RETURN IF(ITER.EQ.ITMAX)THEN WRITE(3,*)' Rotation has not converged !' IERR=1 RETURN ENDIF 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,XS,XB,NAT) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR,A,XS,XB,NAT) 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,XS,XB,NAT) 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,XS,XB,NAT) 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,XS,XB,NAT) IMPLICIT none integer*4 NAT,I real*8 A(3,3),XS(NAT,3),XB(NAT,3),S1,S2,S3,C1,C2,C3,R,FU, 1TX,TY,TZ,ANG(*) 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 R=0.0d0 DO 1 I=1,NAT 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 R=R+(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=R RETURN END subroutine center(X,n) implicit none integer*4 n,ia,ix real*8 X(n,3),t do 1 ix=1,3 t=0.0d0 do 2 ia=1,n 2 t=t+X(ia,ix) t=t/dble(n) do 1 ia=1,n 1 X(ia,ix)=X(ia,ix)-t return end SUBROUTINE READS(f,N3,S,E,NM) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION S(N3,N3),E(*) character*(*) f SFACR=0.02342179d0 OPEN(4,FILE=f,STATUS='OLD') read(4,*)NM,nat,nat do 1 i=1,NAT 1 read(4,*) read(4,*) DO 2 I=1,NAT DO 2 J=1,NM 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,NM) 4000 FORMAT(6F11.6) CLOSE(4) do 4 i=1,nm do 4 k=1,n3 4 S(k,i)=S(k,i)*SFACR write(6,*)NM,' modes in '//f RETURN end C SUBROUTINE READFF(f,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) character*(*) f open(20,file=f,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) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) close(20) WRITE(*,*)f//' read' RETURN END c ============================================================ 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 c ============================================================ SUBROUTINE READTEN(s,N,P,A) IMPLICIT none INTEGER*4 N,L,J,I REAL*8 P(N,3,3),A(N,3,3) character*(*) s open(15,file=s) READ (15,*) DO 10 L=1,N DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) DO 220 L=1,N DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) close(15) WRITE(6,*)s//' read' RETURN END c ============================================================ SUBROUTINE PO(NAT,P,A,AR,X) C Put axial tensor to position X implicit none integer*4 NAT,I,L,J,ID,IG real*8 P(NAT,3,3),A(NAT,3,3),R(3),AR(NAT,3,3),X(*),BOHR,EPS BOHR=0.52917705993d0 DO 3 I=1,3 3 R(I)=X(I)/BOHR DO 1 L=1,NAT DO 1 J=1,3 DO 1 I=1,3 AR(L,J,I)=A(L,J,I) DO 1 ID=1,3 DO 1 IG=1,3 1 AR(L,J,I)=AR(L,J,I)-0.25d0*EPS(I,IG,ID)*R(IG)*P(L,J,ID) WRITE(6,6)X(1),X(2),X(3) 6 format(' AAT shifted to ',3F12.2) RETURN END c ============================================================ subroutine dd(u,m,P,AAT,NAT,nm,E,im,DS,RS,SB,io,j) implicit none integer*4 NAT,ix,ia,iy,nm,im,io,j real*8 u(3),m(3),P(NAT,3,3),AAT(NAT,3,3),SB(3*NAT,nm),E(*), 1FD,FC,DS,RS,RR,c,pi,SFACR pi=4.0d0*atan(1.0d0) FD =SQRT(388.9219d0) FC=131.14D-8 SFACR=0.02342179d0 do 13 ix=1,3 u(ix)=0.0d0 m(ix)=0.0d0 do 131 ia=1,NAT do 131 iy=1,3 u(ix)=u(ix)+SB(iy+3*(ia-1),j)* P(ia,iy,ix) 131 m(ix)=m(ix)+SB(iy+3*(ia-1),j)*AAT(ia,iy,ix) u(ix)=u(ix)/dsqrt(E(im))*FD/SFACR 13 m(ix)=m(ix)*dsqrt(E(im))*FC/SFACR DS=u(1)*u(1)+u(2)*u(2)+u(3)*u(3) RS=u(1)*m(1)+u(2)*m(2)+u(3)*m(3) RR=m(1)*m(1)+m(2)*m(2)+m(3)*m(3) write(6,6002)io,im,E(im),DS,RS,u,m 6002 format('over:',i3,' mode:',i3,f10.2,' D:',g10.4,' R:',G10.4, 1'U,M:',2(3G10.4,' ')) c magnetic dipole part projected to u: c=RS/dsqrt(RR*DS) write(6,6003)c*100.0d0,acos(c)*180.0d0/pi 6003 format(f10.2,'% (',f10.2,' deg)') return end