c rmsdnmrv2.f
c
c generate rmsd and order parameter 
c    from md generated pdb ensemble files
c This program is for a single chain (A chain)
c Ignore all B chain lines
c
c michael h. peters september 25, 2023
c
c *.pdb files are created using VMD. Do not save first 50 snapshots.
c Change HSD to HIS. Change last entry to PROA
c Delete water and ion atom lines to save space 
c Use for example in unix:
c grep -v "TIP3S" oldfile.pdb > newfile.pdb
c or Notepad++
c 
c DO NOT delete END statements since these are used as a string
c for ensemble identification
c
       integer nes
c
c maximum number of snapshots from file
       nes=250
c
       call disp(nes)
       end
c
       subroutine disp(nes)
c
c 1200 total residues max and 300 hundred ensembles
       implicit real(a-h,o-z)
       real xcaa(1200,300),xcab(1200,300),xcavga(1200),xcavgb(1200)
       real ycaa(1200,300),ycab(1200,300),ycavga(1200),ycavgb(1200)
       real zcaa(1200,300),zcab(1200,300),zcavga(1200),zcavgb(1200)
       real rmsda(1200),rmsdb(1200),ca(1200,300)
       real xna(1200,300),yna(1200,300),zna(1200,300)
       real xca(1200,300),yca(1200,300),zca(1200,300)
       real xha(1200,300),yha(1200,300),zha(1200,300)
       real xnha(1200,300),ynha(1200,300),znha(1200,300),rnha(1200,300)
       real xnhpa(1200,300),ynhpa(1200,300),znhpa(1200,300)
       real rnhpa(1200,300)
       real xnb(1200,300),ynb(1200,300),znb(1200,300)
       real xhb(1200,300),yhb(1200,300),zhb(1200,300)
       real xnhb(1200,300),ynhb(1200,300),znhb(1200,300),rnhb(1200,300)
       real xnhpb(1200,300),ynhpb(1200,300),znhpb(1200,300),cb(1200,300)
c
       real xcabp(1200,300),ycabp(1200,300),zcabp(1200,300)
c
       real xcaao(300),ycaao(300),zcaao(300)
       real xhaao(300),yhaao(300),zhaao(300)
       real xcao(300),ycao(300),zcao(300)
       real xnao(300),ynao(300),znao(300)
       real xhaa(1200,300),yhaa(1200,300),zhaa(1200,300)
c
       real x1x2a(300),y1y2a(300),z1z2a(300),x1x2b(300),y1y2b(300)
       real z1z2b(300),thetaa(300),thetab(300)
c
       real xsum(20000),ysum(20000),zsum(20000)
       real A(3,3),CM(3,3),RT(3,3),RTP(3,3),B(3,3)
       real RTPI(3,3),CMT(3,3)
c
       integer resnr,nbar,svrsna(1200),svrsnb(1200)
       integer ies,iesmax,nes,k
c
       integer na,nb,nseqa,nseqb
       character*6 card
       character*1 alt,chid,code
       character*1 chida,chidb
       character*4 atomna,proi
       character*3 resna
c
       open(unit=17, file='7thkOutm50.pdb', status='old')
       open(unit=20, file='ca.pdb', status='new')
       open(unit=21, file='rot.pdb', status='new')
       open(unit=31, file='rotri.pdb', status='new')
       open(unit=32, file='avg.pdb', status='new')
       open(unit=33, file='trans.pdb', status='new')
       open(unit=34, file='dummy.pdb', status='new')
c
c put in new chain id identification
c ignore B chain lines or remove them altogether
c
c
       chida="A"
       chidb="B"
c
c identify c-alpha, n, and h-n atom positions
c note that iesmax is the maximum time interval
c for time correlation and must be less than nes
c its value must be studied for each system to ensure accuracy
c
c fort.16 is rmsfa fort.18 is rmsfb
c fort.22 is rota fort.26 is rotb
c fort.24 is s2a fort.28 is s2b
c
       rewind(17)
       ies=1
       iesmax=150
       nseqa=0
       nseqb=0     
c
c identify body frame origin 
c
       iout=1
       do
       read(unit=17,fmt=12,end=98,iostat=ios)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge,proi
 12    format(a6,i5,1x,a4,a1,a3,1x,a1,i4,a1,3x,3f8.3,2x,f4.2,2x,
     1             f4.2,6x,a4)
       if(resnr.eq.398.and.atomna.eq." CA ") then
       xcaao(iout)=x
       ycaao(iout)=y
       zcaao(iout)=z
       end if
       if(resnr.eq.398.and.atomna.eq." C  ") then
       xcao(iout)=x
       ycao(iout)=y
       zcao(iout)=z
       end if
       if(resnr.eq.398.and.atomna.eq." N  ") then
       xnao(iout)=x
       ynao(iout)=y
       znao(iout)=z
       end if
       if(resnr.eq.398.and.atomna.eq." HA ") then
       xhaao(iout)=x
       yhaao(iout)=y
       zhaao(iout)=z
       end if
c
       if(card.eq."END   ") iout=iout+1
c
       end do
c
 98    continue
c
       rewind(17)
       rewind(33)
c
       write(unit=33,fmt=801)ies
 801   format('MODEL ',4x,i4)
c
       do 
c
       read(unit=17,fmt=12,end=99,iostat=ios)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge,proi
c 12    format(a6,i5,1x,a4,a1,a3,1x,a1,i4,a1,3x,3f8.3,2x,f4.2,2x,
c     1             f4.2,6x,a4)
c
c write to trans.pdb
c
       if(card.eq."END   ")then
       write(unit=33,fmt=800)
 800   format('ENDMDL')
       iesp1=ies+1
       write(unit=33,fmt=801)iesp1
       else
       xp=x-xcaao(ies)
       yp=y-ycaao(ies)
       zp=z-zcaao(ies)
       write(unit=33,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,xp,yp,zp,amass,acharge,proi
       end if
c
c average last ten snapshot locations relative to body origin avg.pdb
       if(ies.eq.(nes-10))then
       xsum(nbar)=0.0
       ysum(nbar)=0.0
       zsum(nbar)=0.0
       end if
c
       if(ies.ge.(nes-9).and.ies.le.nes)then
       xsum(nbar)=xsum(nbar)+(x-xcaao(ies))
       ysum(nbar)=ysum(nbar)+(y-ycaao(ies))
       zsum(nbar)=zsum(nbar)+(z-zcaao(ies))
       end if
c
       if(ies.eq.nes)then
       xavg=xsum(nbar)/10.0
       yavg=ysum(nbar)/10.0
       zavg=zsum(nbar)/10.0
       write(unit=32,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,xavg,yavg,zavg,amass,acharge,proi
       end if
c
c compute ca.pdb body frame for illustration(last 50 to go along with avg.pdb)
c       if(resna.eq."PRO")go to 94
c
       if(proi.eq."PROA".and.atomna.eq." CA ")then
c nseqa numbers the residues from 1 to the max number
       nseqa=nseqa+1
       xcaa(resnr,ies)=x-xcaao(ies)
       ycaa(resnr,ies)=y-ycaao(ies)
       zcaa(resnr,ies)=z-zcaao(ies)
       svrsna(nseqa)=resnr
       chid=chida
c       write(23,200)resnr,ies,xcaa(resnr,ies)
c 200   format('resnr',i6,'ies',i3,'xcaa',e15.8)
c
       if(ies.gt.200)then
       write(unit=20,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,xcaa(resnr,ies),ycaa(resnr,ies),zcaa(resnr,ies),
     1     amass,acharge
       end if
       end if
c
c now compute rot.pdb body frame for illustration
      if(proi.eq."PROA".and.(atomna.eq." HN ".or.atomna.eq." H  "
     1                           .or.atomna.eq." HT1"))then
       xha(resnr,ies)=x-xcaao(ies)
       yha(resnr,ies)=y-ycaao(ies)
       zha(resnr,ies)=z-zcaao(ies)
       chid=chida
       if(ies.gt.200)then
       write(unit=21,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,xha(resnr,ies),yha(resnr,ies),zha(resnr,ies),
     1     amass,acharge
       end if
       end if
c
       if(proi.eq."PROA".and.atomna.eq." N  ")then
       xna(resnr,ies)=x-xcaao(ies)
       yna(resnr,ies)=y-ycaao(ies)
       zna(resnr,ies)=z-zcaao(ies)
       chid=chida
       if(ies.gt.200)then
       write(unit=21,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,xna(resnr,ies),yna(resnr,ies),zna(resnr,ies),
     1     amass,acharge
       end if
       end if
c
       if(proi.eq."PROA".and.atomna.eq." C  ")then
       xca(resnr,ies)=x-xcaao(ies)
       yca(resnr,ies)=y-ycaao(ies)
       zca(resnr,ies)=z-zcaao(ies)
       chid=chida
       end if
c
       if(proi.eq."PROA".and.atomna.eq." HA ")then
       xhaa(resnr,ies)=x-xcaao(ies)
       yhaa(resnr,ies)=y-ycaao(ies)
       zhaa(resnr,ies)=z-zcaao(ies)
       chid=chida
       end if
c
       if(proi.eq."PROB".and.atomna.eq." CA ")then
c nseqb numbers the residues from 1 to the max number
       nseqb=nseqb+1
       xcab(resnr,ies)=x
       ycab(resnr,ies)=y
       zcab(resnr,ies)=z
       svrsnb(nseqb)=resnr
       chid=chidb
c       write(23,200)resnr,ies,xcab(resnr,ies)
c
      write(unit=20,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge
      end if
c
      if(proi.eq."PROB".and.(atomna.eq." HN ".or.atomna.eq." H  "
     1                           .or.atomna.eq." HT1"))then
       xhb(resnr,ies)=x
       yhb(resnr,ies)=y
       zhb(resnr,ies)=z
       chid=chidb
       write(unit=21,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge
       end if
c
       if(proi.eq."PROB".and.atomna.eq." N  ")then
       xnb(resnr,ies)=x
       ynb(resnr,ies)=y
       znb(resnr,ies)=z
       chid=chidb
       write(unit=21,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge
       end if
c
       if(card.eq."END   ")then
       ies=ies+1
       if(ies.gt.nes)go to 99
       nseqa=0
       nseqb=0
       end if
c
 94    continue
c
       end do
c
 99    continue
c
       ies=ies-1
       write(23,201)nseqa,nseqb,ies
 201   format('number of A residues',i6,
     1        'number of B residues',i6,'number of ensembles',i6)
c
c now compute the average position for each c-alpha body frame
c
      do 22 i=1,nseqa
      resnr=svrsna(i)
      xcavga(resnr)=0.0
      ycavga(resnr)=0.0
      zcavga(resnr)=0.0
c
      do 23 j=1,nes
      xcavga(resnr)=xcavga(resnr)+xcaa(resnr,j)
      ycavga(resnr)=ycavga(resnr)+ycaa(resnr,j)
      zcavga(resnr)=zcavga(resnr)+zcaa(resnr,j)
 23   continue
      xcavga(resnr)=xcavga(resnr)/real(nes)
      ycavga(resnr)=ycavga(resnr)/real(nes)
      zcavga(resnr)=zcavga(resnr)/real(nes)
c      write(23,202)resnr,xcavga(resnr),ycavga(resnr),zcavga(resnr)
c 202  format('resnr',i6,'xyzcavga',3e15.8)
c
 22   continue
c
      do 24 i=1,nseqb
      resnr=svrsnb(i)
      xcavgb(resnr)=0.0
      ycavgb(resnr)=0.0
      zcavgb(resnr)=0.0
c
      do 25 j=1,nes
      xcabp(resnr,j)=xcab(resnr,j)-xcab(433,j)
      ycabp(resnr,j)=ycab(resnr,j)-ycab(433,j)
      zcabp(resnr,j)=zcab(resnr,j)-zcab(433,j)
c
      xcavgb(resnr)=xcavgb(resnr)+xcabp(resnr,j)
      ycavgb(resnr)=ycavgb(resnr)+ycabp(resnr,j)
      zcavgb(resnr)=zcavgb(resnr)+zcabp(resnr,j)
 25   continue
      xcavgb(resnr)=xcavgb(resnr)/real(nes)
      ycavgb(resnr)=ycavgb(resnr)/real(nes)
      zcavgb(resnr)=zcavgb(resnr)/real(nes)
c      write(23,203)resnr,xcavgb(resnr),ycavgb(resnr),zcavgb(resnr)
c 203  format('resnr',i6,'xyzcavgb',3e15.8)
c
 24   continue
c
c now compute the rmsd for each c-alpha in body frame
c
      do 26 i=1,nseqa
      resnr=svrsna(i)
      rmsda(resnr)=0.0
c
      do 27 j=1,nes
      rmsda(resnr)=rmsda(resnr)+((xcaa(resnr,j)-xcavga(resnr))**2
     1                          +(ycaa(resnr,j)-ycavga(resnr))**2
     1                          +(zcaa(resnr,j)-zcavga(resnr))**2)
 27   continue
      rmsda(resnr)=sqrt(rmsda(resnr)/real(nes))
      write(16,204)resnr,rmsda(resnr)
 204  format(i6,1x,e15.8)
c
 26   continue
c
      do 28 i=1,nseqb
      resnr=svrsnb(i)
      rmsdb(resnr)=0.0
c
      do 29 j=1,nes
      rmsdb(resnr)=rmsdb(resnr)+((xcabp(resnr,j)-xcavgb(resnr))**2
     1                          +(ycabp(resnr,j)-ycavgb(resnr))**2
     1                          +(zcabp(resnr,j)-zcavgb(resnr))**2)
 29   continue
      rmsdb(resnr)=sqrt(rmsdb(resnr)/real(nes))
      write(18,204)resnr,rmsdb(resnr)
c
 28   continue
c
c correct pdb file for tumbling
c initiate rotri.pdb
       rewind(33)
       rewind(31)
       do
       read(unit=33,fmt=12,end=100,iostat=ios)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge,proi
       write(unit=31,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge,proi
       end do
 100   continue
c
       nesm1=nes-1
       do 66 ies=1,nesm1
       rewind(31)
       rewind(34)
c determine rotational transformation matrix for each ies
       RT(1,1)=xca(398,ies)
       RT(1,2)=yca(398,ies)
       RT(1,3)=zca(398,ies)
       RT(2,1)=xhaa(398,ies)
       RT(2,2)=yhaa(398,ies)
       RT(2,3)=zhaa(398,ies)
       RT(3,1)=xna(398,ies)
       RT(3,2)=yna(398,ies)
       RT(3,3)=zna(398,ies)
c 
c       write(44,*)RT
       RTP(1,1)=xca(398,ies+1)
       RTP(1,2)=yca(398,ies+1)
       RTP(1,3)=zca(398,ies+1)
       RTP(2,1)=xhaa(398,ies+1)
       RTP(2,2)=yhaa(398,ies+1)
       RTP(2,3)=zhaa(398,ies+1)
       RTP(3,1)=xna(398,ies+1)
       RTP(3,2)=yna(398,ies+1)
       RTP(3,3)=zna(398,ies+1)
c
       A=RTP
c       write(44,*)A
c Calculate the inverse determinant of the matrix
       detinv = 1.0/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)
     1         - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)
     1         + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1))

c Calculate the inverse of the matrix
       B(1,1) = +detinv * (A(2,2)*A(3,3) - A(2,3)*A(3,2))
       B(2,1) = -detinv * (A(2,1)*A(3,3) - A(2,3)*A(3,1))
       B(3,1) = +detinv * (A(2,1)*A(3,2) - A(2,2)*A(3,1))
       B(1,2) = -detinv * (A(1,2)*A(3,3) - A(1,3)*A(3,2))
       B(2,2) = +detinv * (A(1,1)*A(3,3) - A(1,3)*A(3,1))
       B(3,2) = -detinv * (A(1,1)*A(3,2) - A(1,2)*A(3,1))
       B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2))
       B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1))
       B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1))
c
       RTPI=B
c       write(44,*)B
c
       CMT=MATMUL(RTPI,RT)
c
       CM=TRANSPOSE(CMT)
       write(44,*)CM
c
       j=1
       do 
       read(unit=31,fmt=12,end=101,iostat=ios)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge,proi
c
       if(j.le.ies)then
       xnew=x
       ynew=y
       znew=z
       else
       xnew=CM(1,1)*x+CM(1,2)*y+CM(1,3)*z
       ynew=CM(2,1)*x+CM(2,2)*y+CM(2,3)*z
       znew=CM(3,1)*x+CM(3,2)*y+CM(3,3)*z
       end if
c
c write to dummy pdb file
c
      if(card.eq."ENDMDL")then
      j=j+1
      write(unit=34,fmt=800)
      else
      write(unit=34,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,xnew,ynew,znew,amass,acharge,proi
      end if
      end do
c
 101  continue
c
c now update 31
      rewind(34) 
      rewind(31)
      iout=1
      do
      read(unit=34,fmt=12,end=102,iostat=ios)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge,proi
      write(unit=31,fmt=12)card,nbar,atomna,alt,
     1     resna,chid,resnr,
     1     code,x,y,z,amass,acharge,proi
c update atom positions
      if(proi.eq."PROA".and.(atomna.eq." HN ".or.atomna.eq." H  "
     1                           .or.atomna.eq." HT1"))then
       xha(resnr,iout)=x
       yha(resnr,iout)=y
       zha(resnr,iout)=z
       end if
c
       if(proi.eq."PROA".and.atomna.eq." N  ")then
       xna(resnr,iout)=x
       yna(resnr,iout)=y
       zna(resnr,iout)=z
       end if
c
       if(proi.eq."PROA".and.atomna.eq." C  ")then
       xca(resnr,iout)=x
       yca(resnr,iout)=y
       zca(resnr,iout)=z
       end if
c
      if(proi.eq."PROA".and.(atomna.eq." HA "))then
       xhaa(resnr,iout)=x
       yhaa(resnr,iout)=y
       zhaa(resnr,iout)=z
       end if
c
      if(card.eq."ENDMDL") iout=iout+1
      end do
 102  continue
c
 66   continue
c
c       iout=1
c       rewind(31)
c       do
c       read(unit=31,fmt=12,end=97,iostat=ios)card,nbar,atomna,alt,
c     1     resna,chid,resnr,
c     1     code,x,y,z,amass,acharge,proi
c
c       if(proi.eq."PROA".and.atomna.eq." N  ") then
c       xna(resnr,iout)=x
c       yna(resnr,iout)=y
c       zna(resnr,iout)=z
c       end if
c       if(proi.eq."PROA".and.atomna.eq." HN ".or.atomna.eq." H  "
c     1                           .or.atomna.eq." HT1") then
c       xha(resnr,iout)=x
c       yha(resnr,iout)=y
c       zha(resnr,iout)=z
c       end if
c
c       if(card.eq."ENDMDL") iout=iout+1
c       end do
c
c 97    continue
c
c compute unit vectors for order parameters
c
      do 38 i=1,nseqa
      resnr=svrsna(i)
c
      do 39 j=1,nes
      xnha(resnr,j)=xha(resnr,j)-xna(resnr,j)
      ynha(resnr,j)=yha(resnr,j)-yna(resnr,j)
      znha(resnr,j)=zha(resnr,j)-zna(resnr,j)
      rnha(resnr,j)=sqrt(xnha(resnr,j)**2+ynha(resnr,j)**2
     1                     +znha(resnr,j)**2)
c
      xnhpa(resnr,j)=xnha(resnr,j)/rnha(resnr,j)
      ynhpa(resnr,j)=ynha(resnr,j)/rnha(resnr,j)
      znhpa(resnr,j)=znha(resnr,j)/rnha(resnr,j)
c
 39   continue
c
 38   continue
c
      do 40 i=1,nseqb
      resnr=svrsnb(i)
c
      do 41 j=1,nes
      xnhb(resnr,j)=xhb(resnr,j)-xnb(resnr,j)
      ynhb(resnr,j)=yhb(resnr,j)-ynb(resnr,j)
      znhb(resnr,j)=zhb(resnr,j)-znb(resnr,j)
      rnhb(resnr,j)=sqrt(xnhb(resnr,j)**2+ynhb(resnr,j)**2
     1                     +znhb(resnr,j)**2)
c
      xnhpb(resnr,j)=xnhb(resnr,j)/rnhb(resnr,j)
      ynhpb(resnr,j)=ynhb(resnr,j)/rnhb(resnr,j)
      znhpb(resnr,j)=znhb(resnr,j)/rnhb(resnr,j)
 41   continue
c
 40   continue
c
c for each residue determine its order parameter
c iesmax is the maximum ensemble number for averaging 
c which must be less that the total number of time snapshots nes
c
      do 50 i=1,nseqa
      resnr=svrsna(i)
c
      do 51 j=1,iesmax
      ca(resnr,j)=0.0
      k=nes-j
      realk=real(k)
      do 52 m=1,k
c      write(18,207)i,j,k,m,ca(i,j)
c 207  format('i j k m',4i5,e15.8)
      mpj=m+j
      ca(resnr,j)=ca(resnr,j)+(1.0/2.0)*(3.0*
     1           ((xnhpa(resnr,m)*xnhpa(resnr,mpj)
     1            +ynhpa(resnr,m)*ynhpa(resnr,mpj)
     1            +znhpa(resnr,m)*znhpa(resnr,mpj))**2.0)-1.0)/realk
 52   continue
c
      write(22,206)resnr,j,ca(resnr,j)
 206  format(2i5,2x,e15.8)
c next interval(j) 
 51   continue
c next residue
c
c write order parameter as last value
      write(24,208)resnr,ca(resnr,iesmax)
 208  format(i5,2x,e15.8)
c
 50   continue
c
      do 53 i=1,nseqb
      resnr=svrsnb(i)
c
      do 54 j=1,iesmax
      cb(resnr,j)=0.0
      k=nes-j
      realk=real(k)
      do 55 m=1,k
c      write(18,207)i,j,k,m,ca(i,j)
c 207  format('i j k m',4i5,e15.8)
      mpj=m+j
      cb(resnr,j)=cb(resnr,j)+(1.0/2.0)*(3.0*
     1           ((xnhpb(resnr,m)*xnhpb(resnr,mpj)
     1            +ynhpb(resnr,m)*ynhpb(resnr,mpj)
     1            +znhpb(resnr,m)*znhpb(resnr,mpj))**2.0)-1.0)/realk
 55   continue
c
      write(26,206)resnr,j,cb(resnr,j)
c next interval(j) 
 54   continue
c next residue
c
c write order parameter as last value
      write(28,208)resnr,cb(resnr,iesmax)
c
 53   continue
c
 988  continue
c
      close(16)
      close(17)
      close(18)
      close(20)
      close(21)
      close(22)
      close(23)
      close(24)
      close(26)
      close(28)  
      close(31)
      close(32)  
      close(33)
      close(34)        
c  
      end subroutine disp
     









