program nmaai implicit none integer*4 NAT,ia,io,ix,nao,NATS,IERR,N3,NM,N,i,j,is,jx,m,i1,i2 integer*4,allocatable::ND(:),IBONDSB(:,:),NBOB(:),INDBIG(:,:), 1NDS(:),IBONDS(:,:),NBO(:) real*8,allocatable::R(:,:),XS(:,:),XB(:,:),RS(:,:),S(:,:),E(:), 1F(:,:),SB(:,:),EB(:) real*8 A(3,3),sx,sy,sz,fi(2,2),w1,w2,CM character*8 s8 write(6,6009) 6009 format(/,' Gets interaction for AI modes based on NMA ',/, 1' in a peptide',/, 2' Input: OVER.TXT',/, 3' BIG.X',/, 5' BIG.FC',/, 4' SMALL.X',/, 4' SMALLF.INP',/) if(iargc().gt.0)then call getarg(1,s8) read(s8,*)m else m=8 endif write(s8,888)m 888 format(i8) write(6,*)'mode ',m,' selected' 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(2,NAT)) call RDOVER(INDBIG,NAT,NDS,ND) N3=3*NATS allocate(S(N3,N3),E(N3)) call READS('SMALLF.INP',N3,S,E,NM) N=3*NAT allocate(F(N,N)) call READFF('BIG.FC',N,F) allocate(SB(N,2),EB(2)) do 5 i=1,2 do 5 j=1,N 5 SB(j,i)=0.0d0 do 1 io=1,2 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 sx=A(1,1)*S(3*is-2,m)+A(1,2)*S(3*is-1,m)+A(1,3)*S(3*is,m) sy=A(2,1)*S(3*is-2,m)+A(2,2)*S(3*is-1,m)+A(2,3)*S(3*is,m) sz=A(3,1)*S(3*is-2,m)+A(3,2)*S(3*is-1,m)+A(3,3)*S(3*is,m) SB(3*ia-2,io)=sx SB(3*ia-1,io)=sy SB(3*ia ,io)=sz endif 4 continue EB(io)=E(m) write(6,*)' A:' 1 write(6,608)A 608 format(3f12.6) call wrf('CONTROL.F.INP',SB,R,NAT,2,EB,ND) do 7 i1=1,len(s8) 7 if(s8(i1:i1).ne.' ')goto 8 8 do 9 i2=len(s8),1,-1 9 if(s8(i2:i2).ne.' ')goto 10 10 continue do 6 i=1,2 do 6 j=1,2 fi(i,j)=0.0d0 do 6 ix=1,N do 6 jx=1,N 6 fi(i,j)=fi(i,j)+SB(ix,i)*F(ix,jx)*SB(jx,j) CM=219474.0d0 w1=dsqrt(fi(1,1)) w2=dsqrt(fi(2,2)) write(6,607)w1*CM,fi(1,2)/dsqrt(w1*w2)*CM/2.0d0,w2*CM 607 format(' H11: ',F12.2,' H12: ',f12.2,' H22: ',f12.2) open(7,file='v12.'//s8(i1:i2)) write(7,700)fi(1,2)/dsqrt(w1*w2)*CM/2.0d0 700 format(2f12.2) close(7) open(7,file='v11v22.'//s8(i1:i2)) write(7,700)w1*219474.0d0,w2*219474.0d0 close(7) 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) implicit none integer*4 NAT,INDBIG(2,NAT),I,IAB,NDS(*),ND(*) open(9,file='OVER.TXT',status='old') do 1 I=1,2 READ(9,*)(INDBIG(I,IAB),IAB=1,NAT) 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