triss; AFM; Lorentz restrains included -debug might be on
[unres4.git] / source / unres / io_base.f90
index f5c7bbf..0e1a986 100644 (file)
@@ -5,10 +5,18 @@
       implicit none
 !-----------------------------------------------------------------------------
 ! Max. number of AA residues
-      integer,parameter :: maxres=4000!1200
+      integer,parameter :: maxres=6000!1200
 ! Appr. max. number of interaction sites
       integer,parameter :: maxres2=2*maxres
+!      parameter (maxres6=6*maxres)
+!      parameter (mmaxres2=(maxres2*(maxres2+1)/2))
 !-----------------------------------------------------------------------------
+! Max. number of S-S bridges
+!      integer,parameter :: maxss=20
+!-----------------------------------------------------------------------------
+! Max. number of derivatives of virtual-bond and side-chain vectors in theta
+! or phi.
+!      integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2
 !-----------------------------------------------------------------------------
 !
 !
@@ -60,6 +68,7 @@
         write(iout,*) ' iss:',(iss(i),i=1,ns)
 ! Check whether the specified bridging residues are cystines.
       do i=1,ns
+         write(iout,*) i,iss(i)
        if (itype(iss(i)).ne.1) then
          if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
          'Do you REALLY think that the residue ',&
 #endif
         endif
       enddo
+      if (dyn_ss) then 
+        if(.not.allocated(ihpb)) allocate(ihpb(ns))
+        if(.not.allocated(jhpb)) allocate(jhpb(ns))
+      endif
 ! Read preformed bridges.
       if (ns.gt.0) then
-      read (inp,*) nss
+        read (inp,*) nss
       if (nss.gt.0) then
         if(.not.allocated(ihpb)) allocate(ihpb(nss))
         if(.not.allocated(jhpb)) allocate(jhpb(nss))
         enddo
       endif
       endif
-!write(iout,*) "end read_bridge"
+!      write(iout,*) "end read_bridge"
       return
       end subroutine read_bridge
 !-----------------------------------------------------------------------------
 !      include 'COMMON.CONTROL'
 !      include 'COMMON.LOCAL'
 !      include 'COMMON.INTERACT'
+! Read coordinates from input
 !
 !el local variables
       integer :: l,k,j,i,kanal
 
-! Read coordinates from input
-!
       read(kanal,'(8f10.5)',end=10,err=10) &
         ((c(l,k),l=1,3),k=1,nres),&
         ((c(l,k+nres),l=1,3),k=nnt,nct)
       end subroutine read_threadbase
 !-----------------------------------------------------------------------------
 #ifdef WHAM_RUN
-      subroutine read_angles(kanal,iscor,energ,iprot,*)
+!el      subroutine read_angles(kanal,iscor,energ,iprot,*)
+      subroutine read_angles(kanal,*)
 
       use geometry_data
       use energy_data
       subroutine read_angles(kanal,*)
 
       use geometry_data
+   !  use energy
+   !  use control
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.GEO'
       card=card(:ilen(card)+1)//karta
       return
       end subroutine card_concat
+!----------------------------------------------------------------------------
+      subroutine read_afminp
+      use geometry_data
+      use energy_data
+      use control_data, only:out1file
+      use MPI_data
+      character*320 afmcard
+      integer i
+      print *, "wchodze"
+      call card_concat(afmcard,.true.)
+      call readi(afmcard,"BEG",afmbeg,0)
+      call readi(afmcard,"END",afmend,0)
+      call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
+      call reada(afmcard,"VEL",velAFMconst,0.0d0)
+      print *,'FORCE=' ,forceAFMconst
+!------ NOW PROPERTIES FOR AFM
+       distafminit=0.0d0
+       do i=1,3
+        distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
+       enddo
+        distafminit=dsqrt(distafminit)
+        print *,'initdist',distafminit
+      return
+      end subroutine read_afminp
 !-----------------------------------------------------------------------------
       subroutine read_dist_constr
       use MPI_data
+   !  use control
       use geometry, only: dist
       use geometry_data
       use control_data
 !      write (iout,*) "Calling read_dist_constr"
 !      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 !      call flush(iout)
+      if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
+      if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
+      if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
+      if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
+      if(.not.allocated(forcon)) allocate(forcon(maxdim))
+      if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
+      if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
+      call gen_dist_constr2
+      go to 1712
+      endif
       call card_concat(controlcard,.true.)
       call readi(controlcard,"NFRAG",nfrag_,0)
       call readi(controlcard,"NPAIR",npair_,0)
 !      do i=1,npair_
 !        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 !      enddo
-      if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
-      if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
-      if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
-      if(.not.allocated(forcon)) allocate(forcon(maxdim))
-
+!      if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
+!      if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
+!      if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
+!      if(.not.allocated(forcon)) allocate(forcon(maxdim))
+      
       call flush(iout)
       do i=1,nfrag_
         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
         endif
       enddo 
       do i=1,ndist_
-        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+!        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+!        if (forcon(nhpb+1).gt.0.0d0) then
+!          nhpb=nhpb+1
+!          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+        if (constr_dist.eq.11) then
+        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
+          ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+        else
+!C        print *,"in else"
+        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
+          ibecarb(i),forcon(nhpb+1)
+        endif
         if (forcon(nhpb+1).gt.0.0d0) then
           nhpb=nhpb+1
-          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+          if (ibecarb(i).gt.0) then
+            ihpb(i)=ihpb(i)+nres
+            jhpb(i)=jhpb(i)+nres
+          endif
+          if (dhpb(nhpb).eq.0.0d0) &
+            dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+        endif
+
 #ifdef MPI
           if (.not.out1file .or. me.eq.king) &
           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
            nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
-        endif
       enddo
+ 1712 continue
       call flush(iout)
       return
       end subroutine read_dist_constr
 !-----------------------------------------------------------------------------
+      subroutine gen_dist_constr2
+      use MPI_data
+   !  use control
+      use geometry, only: dist
+      use geometry_data
+      use control_data
+      use energy_data
+      integer :: i,j
+      real(kind=8) :: distance
+      if (constr_dist.eq.11) then
+             do i=nstart_sup,nstart_sup+nsup-1
+              do j=i+2,nstart_sup+nsup-1
+                 distance=dist(i,j)
+                 if (distance.le.15.0) then
+                 nhpb=nhpb+1
+                 ihpb(nhpb)=i+nstart_seq-nstart_sup
+                 jhpb(nhpb)=j+nstart_seq-nstart_sup
+                 forcon(nhpb)=sqrt(0.04*distance)
+                 fordepth(nhpb)=sqrt(40.0/distance)
+                 dhpb(nhpb)=distance-0.1d0
+                 dhpb1(nhpb)=distance+0.1d0
+
+#ifdef MPI
+          if (.not.out1file .or. me.eq.king) &
+          write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
+          nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+          write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
+          nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+            endif
+             enddo
+           enddo
+      else
+      do i=nstart_sup,nstart_sup+nsup-1
+        do j=i+2,nstart_sup+nsup-1
+          nhpb=nhpb+1
+          ihpb(nhpb)=i+nstart_seq-nstart_sup
+          jhpb(nhpb)=j+nstart_seq-nstart_sup
+          forcon(nhpb)=weidis
+          dhpb(nhpb)=dist(i,j)
+        enddo
+      enddo
+      endif
+      return
+      end subroutine gen_dist_constr2
+
+!-----------------------------------------------------------------------------
 #ifdef WINIFL
       subroutine flush(iu)
       return
 
       use geometry_data, only: c,nres
       use energy_data
+   !  use control
       use compare_data
       use MD_data
 !      implicit real*8 (a-h,o-z)
       character*(*) :: tytul
       character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
       integer,dimension(nres) :: ica   !(maxres)
-
+       integer iti1
 !el  local variables
       integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
       real(kind=8) :: etot
       ires=0
       do i=nnt,nct
         iti=itype(i)
+        iti1=itype(i+1)
+        if ((iti.eq.ntyp1).and.(iti1.eq.ntyp1)) cycle
         if (iti.eq.ntyp1) then
           ichain=ichain+1
-          ires=0
+!          ires=0
           write (iunit,'(a)') 'TER'
         else
         ires=ires+1
 ! format.
       use geometry_data, only: c
       use energy_data
+   !  use control
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.CHAIN'
       return
       end subroutine intout
 !-----------------------------------------------------------------------------
+#ifdef CLUSTER
+      subroutine briefout(it,ener,free)!,plik)
+#else
       subroutine briefout(it,ener)
-
+#endif
       use geometry_data
       use energy_data
+   !  use control
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.IOUNITS'
 !     print '(a,i5)',intname,igeom
 !el  local variables
       integer :: i,it
-      real(kind=8) :: ener
-#ifdef WHAM_RUN
-      integer :: iii
-#endif
+      real(kind=8) :: ener,free
+!     character(len=80) :: plik
+!      integer :: iii
 
 #if defined(AIX) || defined(PGI)
       open (igeom,file=intname,position='append')
       open (igeom,file=intname,access='append')
 #endif
 #ifdef WHAM_RUN
-      iii=igeom
+!      iii=igeom
       igeom=iout
 #endif
+      print *,nss
       IF (NSS.LE.9) THEN
+#ifdef CLUSTER
+        WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
+      ELSE
+        WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
+#else
         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)
+#endif
         WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
       ENDIF
 !     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
       end subroutine reads
 #endif
 !-----------------------------------------------------------------------------
+! permut.F
+!-----------------------------------------------------------------------------
+      subroutine permut(isym)
+
+      use geometry_data, only: tabperm
+!      implicit real*8 (a-h,o-z) 
+!      include 'DIMENSIONS'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.VAR'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.GEO'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.CONTROL'
+
+      integer :: n,isym
+!      logical nextp
+!el      external nextp
+      integer,dimension(isym) :: a
+!      parameter(n=symetr)
+!el local variables
+      integer :: kkk,i
+
+      n=isym
+      if (n.eq.1) then
+        tabperm(1,1)=1
+        return
+      endif
+      kkk=0
+      do i=1,n
+      a(i)=i
+      enddo
+   10 print *,(a(i),i=1,n)
+      kkk=kkk+1
+      do i=1,n
+      tabperm(kkk,i)=a(i)
+!      write (iout,*) "tututu", kkk
+      enddo
+      if(nextp(n,a)) go to 10
+      return
+      end subroutine permut
+!-----------------------------------------------------------------------------
+      logical function nextp(n,a)
+
+      integer :: n,i,j,k,t
+!      logical :: nextp
+      integer,dimension(n) :: a
+      i=n-1
+   10 if(a(i).lt.a(i+1)) go to 20
+      i=i-1
+      if(i.eq.0) go to 20
+      go to 10
+   20 j=i+1
+      k=n
+   30 t=a(j)
+      a(j)=a(k)
+      a(k)=t
+      j=j+1
+      k=k-1
+      if(j.lt.k) go to 30
+      j=i
+      if(j.ne.0) go to 40
+      nextp=.false.
+      return
+   40 j=j+1
+      if(a(j).lt.a(i)) go to 40
+      t=a(i)
+      a(i)=a(j)
+      a(j)=t
+      nextp=.true.
+      return
+      end function nextp
+!-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
       end module io_base