Merge branch 'lipid' of mmka.chem.univ.gda.pl:unres into lipid
[unres.git] / source / unres / src_MD-M-newcorr / geomout.F
diff --git a/source/unres/src_MD-M-newcorr/geomout.F b/source/unres/src_MD-M-newcorr/geomout.F
new file mode 100644 (file)
index 0000000..e869b4a
--- /dev/null
@@ -0,0 +1,512 @@
+      subroutine pdbout(etot,tytul,iunit)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.HEADER'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.DISTFIT'
+      include 'COMMON.MD'
+      character*50 tytul
+      character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/
+      dimension ica(maxres)
+      write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
+cmodel      write (iunit,'(a5,i6)') 'MODEL',1
+      if (nhfrag.gt.0) then
+       do j=1,nhfrag
+        iti=itype(hfrag(1,j))
+        itj=itype(hfrag(2,j))
+        if (j.lt.10) then
+           write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') 
+     &           'HELIX',j,'H',j,
+     &           restyp(iti),hfrag(1,j)-1,
+     &           restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
+        else
+             write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') 
+     &           'HELIX',j,'H',j,
+     &           restyp(iti),hfrag(1,j)-1,
+     &           restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
+        endif
+       enddo
+      endif
+
+      if (nbfrag.gt.0) then
+
+       do j=1,nbfrag
+
+        iti=itype(bfrag(1,j))
+        itj=itype(bfrag(2,j)-1)
+
+        write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') 
+     &           'SHEET',1,'B',j,2,
+     &           restyp(iti),bfrag(1,j)-1,
+     &           restyp(itj),bfrag(2,j)-2,0
+
+        if (bfrag(3,j).gt.bfrag(4,j)) then
+
+         itk=itype(bfrag(3,j))
+         itl=itype(bfrag(4,j)+1)
+
+         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,
+     &              2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') 
+     &           'SHEET',2,'B',j,2,
+     &           restyp(itl),bfrag(4,j),
+     &           restyp(itk),bfrag(3,j)-1,-1,
+     &           "N",restyp(itk),bfrag(3,j)-1,
+     &           "O",restyp(iti),bfrag(1,j)-1
+
+        else
+
+         itk=itype(bfrag(3,j))
+         itl=itype(bfrag(4,j)-1)
+
+
+        write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,
+     &              2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') 
+     &           'SHEET',2,'B',j,2,
+     &           restyp(itk),bfrag(3,j)-1,
+     &           restyp(itl),bfrag(4,j)-2,1,
+     &           "N",restyp(itk),bfrag(3,j)-1,
+     &           "O",restyp(iti),bfrag(1,j)-1
+
+
+
+        endif
+         
+       enddo
+      endif 
+
+      if (nss.gt.0) then
+        do i=1,nss
+         if (dyn_ss) then
+          write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)')
+     &         'SSBOND',i,'CYS',idssb(i)-nnt+1,
+     &                    'CYS',jdssb(i)-nnt+1
+         else
+          write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') 
+     &         'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,
+     &                    'CYS',jhpb(i)-nnt+1-nres
+          endif
+        enddo
+      endif
+      
+      iatom=0
+      ichain=1
+      ires=0
+      do i=nnt,nct
+        iti=itype(i)
+        if (iti.eq.ntyp1) then
+          ichain=ichain+1
+          ires=0
+          write (iunit,'(a)') 'TER'
+        else
+        ires=ires+1
+        iatom=iatom+1
+        ica(i)=iatom
+        write (iunit,10) iatom,restyp(iti),chainid(ichain),
+     &     ires,(c(j,i),j=1,3),vtot(i)
+        if (iti.ne.10) then
+          iatom=iatom+1
+          write (iunit,20) iatom,restyp(iti),chainid(ichain),
+     &      ires,(c(j,nres+i),j=1,3),
+     &      vtot(i+nres)
+        endif
+        endif
+      enddo
+      write (iunit,'(a)') 'TER'
+      do i=nnt,nct-1
+        if (itype(i).eq.ntyp1) cycle
+        if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
+          write (iunit,30) ica(i),ica(i+1)
+        else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+          write (iunit,30) ica(i),ica(i+1),ica(i)+1
+        else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
+          write (iunit,30) ica(i),ica(i)+1
+        endif
+      enddo
+      if (itype(nct).ne.10) then
+        write (iunit,30) ica(nct),ica(nct)+1
+      endif
+      do i=1,nss
+C        write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
+c trzeba uporzdkowac 
+C        write (iunit,30) ica(ihpb(i))+1,ica(jhpb(i))+1
+       if (dyn_ss) then
+        write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
+       else
+        write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
+       endif
+      enddo
+      write (iunit,'(a6)') 'ENDMDL'     
+  10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+  20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+  30  FORMAT ('CONECT',8I5)
+      return
+      end
+c------------------------------------------------------------------------------
+      subroutine MOL2out(etot,tytul)
+C Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2 
+C format.
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.HEADER'
+      include 'COMMON.SBRIDGE'
+      character*32 tytul,fd
+      character*3 zahl
+      character*6 res_num,pom,ucase
+#ifdef AIX
+      call fdate_(fd)
+#elif (defined CRAY)
+      call date(fd)
+#else
+      call fdate(fd)
+#endif
+      write (imol2,'(a)') '#'
+      write (imol2,'(a)') 
+     & '#         Creating user name:           unres'
+      write (imol2,'(2a)') '#         Creation time:                ',
+     & fd
+      write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
+      write (imol2,'(a)') tytul
+      write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
+      write (imol2,'(a)') 'SMALL'
+      write (imol2,'(a)') 'USER_CHARGES'
+      write (imol2,'(a)') '\@<TRIPOS>ATOM' 
+      do i=nnt,nct
+        write (zahl,'(i3)') i
+        pom=ucase(restyp(itype(i)))
+        res_num = pom(:3)//zahl(2:)
+        write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
+      enddo
+      write (imol2,'(a)') '\@<TRIPOS>BOND'
+      do i=nnt,nct-1
+        write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
+      enddo
+      do i=1,nss
+        write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
+      enddo
+      write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
+      do i=nnt,nct
+        write (zahl,'(i3)') i
+        pom = ucase(restyp(itype(i)))
+        res_num = pom(:3)//zahl(2:)
+        write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
+      enddo
+  10  FORMAT (I7,' CA      ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
+  30  FORMAT (I7,1x,A,I14,' RESIDUE',I13,' ****  ****')
+      return
+      end
+c------------------------------------------------------------------------
+      subroutine intout
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      include 'COMMON.GEO'
+      include 'COMMON.TORSION'
+      write (iout,'(/a)') 'Geometry of the virtual chain.'
+      write (iout,'(7a)') '  Res  ','         d','     Theta',
+     & '       Phi','       Dsc','     Alpha','      Omega'
+      do i=1,nres
+       iti=itype(i)
+        write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),
+     &     rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),
+     &     rad2deg*omeg(i)
+      enddo
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine briefout(it,ener)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      include 'COMMON.GEO'
+      include 'COMMON.SBRIDGE'
+c     print '(a,i5)',intname,igeom
+#if defined(AIX) || defined(PGI)
+      open (igeom,file=intname,position='append')
+#else
+      open (igeom,file=intname,access='append')
+#endif
+      IF (NSS.LE.9) THEN
+        WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
+      ELSE
+        WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
+        WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
+      ENDIF
+c     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
+      WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
+      WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
+c     if (nvar.gt.nphi+ntheta) then
+        write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
+        write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
+c     endif
+      close(igeom)
+  180 format (I5,F12.3,I2,9(1X,2I3))
+  190 format (3X,11(1X,2I3))
+  200 format (8F10.4)
+      return
+      end
+#ifdef WINIFL
+      subroutine fdate(fd)
+      character*32 fd
+      write(fd,'(32x)')
+      return
+      end
+#endif
+c----------------------------------------------------------------
+#ifdef NOXDR
+      subroutine cartout(time)
+#else
+      subroutine cartoutx(time)
+#endif
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.HEADER'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.DISTFIT'
+      include 'COMMON.MD'
+      double precision time
+#if defined(AIX) || defined(PGI)
+      open(icart,file=cartname,position="append")
+#else
+      open(icart,file=cartname,access="append")
+#endif
+      write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
+C      write (icart,'(i4,$)')
+      if (dyn_ss) then
+       write (icart,'(i4,$)')
+     &   nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)       
+      else
+       write (icart,'(i4,$)')
+     &   nss,(ihpb(j),jhpb(j),j=1,nss)
+       endif 
+       write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,
+     & (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),
+     & (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
+      write (icart,'(8f10.5)')
+     & ((c(k,j),k=1,3),j=1,nres),
+     & ((c(k,j+nres),k=1,3),j=nnt,nct)
+      close(icart)
+      return
+      end
+c-----------------------------------------------------------------
+#ifndef NOXDR
+      subroutine cartout(time)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+      include 'COMMON.SETUP'
+#else
+      parameter (me=0)
+#endif
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.HEADER'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.DISTFIT'
+      include 'COMMON.MD'
+      double precision time
+      integer iret,itmp
+      real xcoord(3,maxres2+2),prec
+
+#ifdef AIX
+      call xdrfopen_(ixdrf,cartname, "a", iret)
+      call xdrffloat_(ixdrf, real(time), iret)
+      call xdrffloat_(ixdrf, real(potE), iret)
+      call xdrffloat_(ixdrf, real(uconst), iret)
+      call xdrffloat_(ixdrf, real(uconst_back), iret)
+      call xdrffloat_(ixdrf, real(t_bath), iret)
+      call xdrfint_(ixdrf, nss, iret) 
+      do j=1,nss
+        call xdrfint_(ixdrf, ihpb(j), iret)
+        call xdrfint_(ixdrf, jhpb(j), iret)
+      enddo
+      call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
+      do i=1,nfrag
+        call xdrffloat_(ixdrf, real(qfrag(i)), iret)
+      enddo
+      do i=1,npair
+        call xdrffloat_(ixdrf, real(qpair(i)), iret)
+      enddo
+      do i=1,nfrag_back
+        call xdrffloat_(ixdrf, real(utheta(i)), iret)
+        call xdrffloat_(ixdrf, real(ugamma(i)), iret)
+        call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
+      enddo
+#else
+      call xdrfopen(ixdrf,cartname, "a", iret)
+      call xdrffloat(ixdrf, real(time), iret)
+      call xdrffloat(ixdrf, real(potE), iret)
+      call xdrffloat(ixdrf, real(uconst), iret)
+      call xdrffloat(ixdrf, real(uconst_back), iret)
+      call xdrffloat(ixdrf, real(t_bath), iret)
+      call xdrfint(ixdrf, nss, iret) 
+      do j=1,nss
+       if (dyn_ss) then
+        call xdrfint(ixdrf, idssb(j)+nres, iret)
+        call xdrfint(ixdrf, jdssb(j)+nres, iret)
+       else
+        call xdrfint(ixdrf, ihpb(j), iret)
+        call xdrfint(ixdrf, jhpb(j), iret)
+       endif 
+      enddo
+      call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
+      do i=1,nfrag
+        call xdrffloat(ixdrf, real(qfrag(i)), iret)
+      enddo
+      do i=1,npair
+        call xdrffloat(ixdrf, real(qpair(i)), iret)
+      enddo
+      do i=1,nfrag_back
+        call xdrffloat(ixdrf, real(utheta(i)), iret)
+        call xdrffloat(ixdrf, real(ugamma(i)), iret)
+        call xdrffloat(ixdrf, real(uscdiff(i)), iret)
+      enddo
+#endif
+      prec=10000.0
+      do i=1,nres
+       do j=1,3
+        xcoord(j,i)=c(j,i)
+       enddo
+      enddo
+      do i=nnt,nct
+       do j=1,3
+        xcoord(j,nres+i-nnt+1)=c(j,i+nres)
+       enddo
+      enddo
+
+      itmp=nres+nct-nnt+1
+#ifdef AIX
+      call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
+      call xdrfclose_(ixdrf, iret)
+#else
+      call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
+      call xdrfclose(ixdrf, iret)
+#endif
+      return
+      end
+#endif
+c-----------------------------------------------------------------
+      subroutine statout(itime)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.HEADER'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.DISTFIT'
+      include 'COMMON.MD'
+      include 'COMMON.SETUP'
+      integer itime
+      double precision energia(0:n_ene)
+      double precision gyrate
+      external gyrate
+      common /gucio/ cm
+      character*256 line1,line2
+      character*4 format1,format2
+      character*30 format
+#ifdef AIX
+      if(itime.eq.0) then
+       open(istat,file=statname,position="append")
+      endif
+#else
+#ifdef PGI
+      open(istat,file=statname,position="append")
+#else
+      open(istat,file=statname,access="append")
+#endif
+#endif
+       if (refstr) then
+         call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
+          write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)')
+     &          itime,totT,EK,potE,totE,
+     &          rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
+          format1="a133"
+        else
+          write (line1,'(i10,f15.2,7f12.3,i5,$)')
+     &           itime,totT,EK,potE,totE,
+     &           amax,kinetic_T,t_bath,gyrate(),me
+          format1="a114"
+        endif
+        if(usampl.and.totT.gt.eq_time) then
+           write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,
+     &      (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),
+     &      (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
+           write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair
+     &             +21*nfrag_back
+        else
+           format2="a001"
+           line2=' '
+        endif
+        if (print_compon) then
+          write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
+     &                                                     ",20f12.3)"
+          write (istat,format) line1,line2,
+     &      (potEcomp(print_order(i)),i=1,nprint_ene)
+        else
+          write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
+          write (istat,format) line1,line2
+        endif
+#if defined(AIX)
+        call flush(istat)
+#else
+        close(istat)
+#endif
+       return
+      end
+c---------------------------------------------------------------                     
+      double precision function gyrate()
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.INTERACT'
+      include 'COMMON.CHAIN'
+      double precision cen(3),rg
+
+      do j=1,3
+       cen(j)=0.0d0
+      enddo
+
+      do i=nnt,nct
+          do j=1,3
+            cen(j)=cen(j)+c(j,i)
+          enddo
+      enddo
+      do j=1,3
+            cen(j)=cen(j)/dble(nct-nnt+1)
+      enddo
+      rg = 0.0d0
+      do i = nnt, nct
+        do j=1,3
+         rg = rg + (c(j,i)-cen(j))**2 
+        enddo
+      end do
+      gyrate = sqrt(rg/dble(nct-nnt+1))
+      return
+      end
+