The itype(;) extension to itype(:,:)
[unres4.git] / source / unres / io_base.f90
index 0a9dc14..303927a 100644 (file)
         write(iout,*) ' iss:',(iss(i),i=1,ns)
 ! Check whether the specified bridging residues are cystines.
       do i=1,ns
-       if (itype(iss(i)).ne.1) then
+         write(iout,*) i,iss(i)
+       if (itype(iss(i),1).ne.1) then
          if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
          'Do you REALLY think that the residue ',&
-          restyp(itype(iss(i))),i,&
+          restyp(itype(iss(i),1)),i,&
          ' can form a disulfide bridge?!!!'
          write (*,'(2a,i3,a)') &
          'Do you REALLY think that the residue ',&
-          restyp(itype(iss(i))),i,&
+          restyp(itype(iss(i),1)),i,&
          ' can form a disulfide bridge?!!!'
 #ifdef MPI
         call MPI_Finalize(MPI_COMM_WORLD,ierror)
 #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
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
           do j=1,3
             dc(j,i+nres)=c(j,i+nres)-c(j,i)
             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
       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
 !      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
 !model      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))
+        iti=itype(hfrag(1,j),1)
+        itj=itype(hfrag(2,j),1)
         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,&
 
        do j=1,nbfrag
 
-        iti=itype(bfrag(1,j))
-        itj=itype(bfrag(2,j)-1)
+        iti=itype(bfrag(1,j),1)
+        itj=itype(bfrag(2,j)-1,1)
 
         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
                  'SHEET',1,'B',j,2,&
 
         if (bfrag(3,j).gt.bfrag(4,j)) then
 
-         itk=itype(bfrag(3,j))
-         itl=itype(bfrag(4,j)+1)
+         itk=itype(bfrag(3,j),1)
+         itl=itype(bfrag(4,j)+1,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,&
 
         else
 
-         itk=itype(bfrag(3,j))
-         itl=itype(bfrag(4,j)-1)
+         itk=itype(bfrag(3,j),1)
+         itl=itype(bfrag(4,j)-1,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)') &
       ichain=1
       ires=0
       do i=nnt,nct
-        iti=itype(i)
-        iti1=itype(i+1)
+        iti=itype(i,1)
+        iti1=itype(i+1,1)
         if ((iti.eq.ntyp1).and.(iti1.eq.ntyp1)) cycle
         if (iti.eq.ntyp1) then
           ichain=ichain+1
       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
+        if (itype(i,1).eq.ntyp1) cycle
+        if (itype(i,1).eq.10 .and. itype(i+1,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
+        else if (itype(i,1).ne.10 .and. itype(i+1,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
+        else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
           write (iunit,30) ica(i),ica(i)+1
         endif
       enddo
-      if (itype(nct).ne.10) then
+      if (itype(nct,1).ne.10) then
         write (iunit,30) ica(nct),ica(nct)+1
       endif
       do i=1,nss
       write (imol2,'(a)') '\@<TRIPOS>ATOM' 
       do i=nnt,nct
         write (zahl,'(i3)') i
-        pom=ucase(restyp(itype(i)))
+        pom=ucase(restyp(itype(i,1)))
         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>SUBSTRUCTURE'
       do i=nnt,nct
         write (zahl,'(i3)') i
-        pom = ucase(restyp(itype(i)))
+        pom = ucase(restyp(itype(i,1)))
         res_num = pom(:3)//zahl(2:)
         write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
       enddo
       write (iout,'(7a)') '  Res  ','         d','     Theta',&
        '       Phi','       Dsc','     Alpha','      Omega'
       do i=1,nres
-       iti=itype(i)
+       iti=itype(i,1)
         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)
 !      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)