c rmsdnmrv2.f
c
c generate rmsd and order parameter 
c    from md generated pdb ensemble files
c
c michael h. peters january 20, 2019
c
c delete water and ion atom lines to save space 
c use:
c grep -v "TIP3S" oldfile.pdb > newfile.pdb
c 
c DO NOT delete END statements since these are used as a string
c for ensemble identification
c
c note that proline residues are skipped everywhere
c
       integer nes
c
c maximum number of snapshots from file after discard
       nes=150
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)
       real xna(1200,300),yna(1200,300),zna(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),ca(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 xcaap(1200,300),ycaap(1200,300),zcaap(1200,300)
       real xcabp(1200,300),ycabp(1200,300),zcabp(1200,300)
c
       real x1x2a(300),y1y2a(300),z1z2a(300),x1x2b(300),y1y2b(300)
       real z1z2b(300),thetaa(300),thetab(300)
       real rmod1a(300),rmod2a(300),rmod1b(300),rmod2b(300)
       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='WTtrimernamdredo2.pdb', status='old')
       open(unit=20, file='ca.pdb', status='new')
       open(unit=21, file='rot.pdb', status='new')
c
c put in correct chain id
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=70
       nseqa=0
       nseqb=0
c
c throw out first 50 snapshots
c
       iout=1
       do
       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
       if(card.eq."END") iout=iout+1
       if(iout.eq.50)go to 98
       end do
c
 98    continue
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
 12    format(a6,i5,1x,a4,a1,a3,1x,a1,i4,a1,3x,3f8.3,2x,f4.2,2x,
     1             f4.2,6x,a4)
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
       ycaa(resnr,ies)=y
       zcaa(resnr,ies)=z
       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
       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."PROA".and.(atomna.eq." HN ".or.atomna.eq." H  "
     1                           .or.atomna.eq." HT1"))then
       xha(resnr,ies)=x
       yha(resnr,ies)=y
       zha(resnr,ies)=z
       chid=chida
       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."PROA".and.atomna.eq." N  ")then
       xna(resnr,ies)=x
       yna(resnr,ies)=y
       zna(resnr,ies)=z
       chid=chida
       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." 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
       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 
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
      xcaap(resnr,j)=xcaa(resnr,j)-xcaa(991,j)
      ycaap(resnr,j)=ycaa(resnr,j)-ycaa(991,j)
      zcaap(resnr,j)=zcaa(resnr,j)-zcaa(991,j)
c
      xcavga(resnr)=xcavga(resnr)+xcaap(resnr,j)
      ycavga(resnr)=ycavga(resnr)+ycaap(resnr,j)
      zcavga(resnr)=zcavga(resnr)+zcaap(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(991,j)
      ycabp(resnr,j)=ycab(resnr,j)-ycab(991,j)
      zcabp(resnr,j)=zcab(resnr,j)-zcab(991,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
c
      do 26 i=1,nseqa
      resnr=svrsna(i)
      rmsda(resnr)=0.0
c
      do 27 j=1,nes
      rmsda(resnr)=rmsda(resnr)+((xcaap(resnr,j)-xcavga(resnr))**2
     1                          +(ycaap(resnr,j)-ycavga(resnr))**2
     1                          +(zcaap(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 compute the hinge angle for each snapshot
c
      do 60 j=1,nes
      x1x2a(j)=(xcaa(405,j)-xcaa(622,j))*(xcaa(991,j)-xcaa(622,j))
      y1y2a(j)=(ycaa(405,j)-ycaa(622,j))*(ycaa(991,j)-ycaa(622,j))
      z1z2a(j)=(zcaa(405,j)-zcaa(622,j))*(zcaa(991,j)-zcaa(622,j))
      rmod1a(j)=sqrt((xcaa(405,j)-xcaa(622,j))**2
     1          +(ycaa(405,j)-ycaa(622,j))**2
     1          +(zcaa(405,j)-zcaa(622,j))**2)
      rmod2a(j)=sqrt((xcaa(991,j)-xcaa(622,j))**2
     1          +(ycaa(991,j)-ycaa(622,j))**2
     1          +(zcaa(991,j)-zcaa(622,j))**2)
      thetaa(j)=acos((x1x2a(j)+y1y2a(j)+z1z2a(j))/(rmod1a(j)*rmod2a(j)))
 60   continue
c
      do 61 j=1,nes
      x1x2b(j)=(xcab(405,j)-xcab(622,j))*(xcab(991,j)-xcab(622,j))
      y1y2b(j)=(ycab(405,j)-ycab(622,j))*(ycab(991,j)-ycab(622,j))
      z1z2b(j)=(zcab(405,j)-zcab(622,j))*(zcab(991,j)-zcab(622,j))
      rmod1b(j)=sqrt((xcab(405,j)-xcab(622,j))**2
     1          +(ycab(405,j)-ycab(622,j))**2
     1          +(zcab(405,j)-zcab(622,j))**2)
      rmod2b(j)=sqrt((xcab(991,j)-xcab(622,j))**2
     1          +(ycab(991,j)-ycab(622,j))**2
     1          +(zcab(991,j)-zcab(622,j))**2)
      thetab(j)=acos((x1x2b(j)+y1y2b(j)+z1z2b(j))/(rmod1b(j)*rmod2b(j)))
c
 61   continue
c
      do 62 j=1,nes
      write(30,350)j,thetaa(j),thetab(j)
 350  format(i5,2f8.3)
      write(31,351)j,rmod1a(j),rmod2a(j),rmod1b(j),rmod2b(j)
 351  format(i5,4e15.8)
 62   continue
c
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)
 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 
      close(16)
      close(17)
      close(18)
      close(20)
      close(21)
      close(22)
      close(23)
      close(24)
      close(26)
      close(28)            
c  
      end subroutine disp