source/unres/io_base.f90

   1       module io_base
   2 !-----------------------------------------------------------------------
   3       use names
   4       use io_units
   5       implicit none
   6 !-----------------------------------------------------------------------------
   7 ! Max. number of AA residues
   8       integer,parameter :: maxres=6000!1200
   9 ! Appr. max. number of interaction sites
  10       integer,parameter :: maxres2=2*maxres
  11 !      parameter (maxres6=6*maxres)
  12 !      parameter (mmaxres2=(maxres2*(maxres2+1)/2))
  13 !-----------------------------------------------------------------------------
  14 ! Max. number of S-S bridges
  15 !      integer,parameter :: maxss=20
  16 !-----------------------------------------------------------------------------
  17 ! Max. number of derivatives of virtual-bond and side-chain vectors in theta
  18 ! or phi.
  19 !      integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2
  20 !-----------------------------------------------------------------------------
  21 !
  22 !
  23 !-----------------------------------------------------------------------------
  24       contains
  25 !-----------------------------------------------------------------------------
  26 ! readrtns_CSA.F
  27 !-----------------------------------------------------------------------------
  28       subroutine read_bridge
  29 ! Read information about disulfide bridges.
  30       use geometry_data, only: nres
  31       use energy_data
  32       use control_data, only:out1file
  33       use MPI_data
  34 !      implicit real*8 (a-h,o-z)
  35 !      include 'DIMENSIONS'
  36 #ifdef MPI
  37       include 'mpif.h'
  38 #endif
  39 !      include 'COMMON.IOUNITS'
  40 !      include 'COMMON.GEO'
  41 !      include 'COMMON.VAR'
  42 !      include 'COMMON.INTERACT'
  43 !      include 'COMMON.LOCAL'
  44 !      include 'COMMON.NAMES'
  45 !      include 'COMMON.CHAIN'
  46 !      include 'COMMON.FFIELD'
  47 !      include 'COMMON.SBRIDGE'
  48 !      include 'COMMON.HEADER'
  49 !      include 'COMMON.CONTROL'
  50 !      include 'COMMON.DBASE'
  51 !      include 'COMMON.THREAD'
  52 !      include 'COMMON.TIME1'
  53 !      include 'COMMON.SETUP'
  54 !el local variables
  55       integer :: i,j,ierror
  56       print *,"ENTER READ"
  57 ! Read bridging residues.
  58       read (inp,*) ns
  59       print *,"ns",ns
  60       if (ns.gt.0) then
  61         allocate(iss(ns))
  62         read (inp,*) (iss(i),i=1,ns)
  63       endif
  64
  65 !      print *,'ns=',ns
  66       if(me.eq.king.or..not.out1file) &
  67         write (iout,*) 'ns=',ns
  68       if (ns.gt.0) &
  69         write(iout,*) ' iss:',(iss(i),i=1,ns)
  70 ! Check whether the specified bridging residues are cystines.
  71       do i=1,ns
  72          write(iout,*) i,iss(i)
  73         if (itype(iss(i),1).ne.1) then
  74           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
  75          'Do you REALLY think that the residue ',&
  76           restyp(itype(iss(i),1),1),i,&
  77          ' can form a disulfide bridge?!!!'
  78           write (*,'(2a,i3,a)') &
  79          'Do you REALLY think that the residue ',&
  80           restyp(itype(iss(i),1),1),i,&
  81          ' can form a disulfide bridge?!!!'
  82 #ifdef MPI
  83          call MPI_Finalize(MPI_COMM_WORLD,ierror)
  84          stop
  85 #endif
  86         endif
  87       enddo
  88       if (dyn_ss) then
  89         if(.not.allocated(ihpb)) allocate(ihpb(ns))
  90         if(.not.allocated(jhpb)) allocate(jhpb(ns))
  91       endif
  92 ! Read preformed bridges.
  93       if (ns.gt.0) then
  94         read (inp,*) nss
  95       if (nss.gt.0) then
  96         if(.not.allocated(ihpb)) allocate(ihpb(nss))
  97         if(.not.allocated(jhpb)) allocate(jhpb(nss))
  98         read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
  99
 100         if(.not.allocated(dhpb)) allocate(dhpb(nss))
 101         if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
 102         if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
 103         if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
 104 ! el Initialize the bridge array
 105         do i=1,nss
 106           dhpb(i)=0.0D0
 107         enddo
 108       endif
 109 !--------------------
 110       if(fg_rank.eq.0) &
 111        write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 112       if (nss.gt.0) then
 113        write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 114         nhpb=nss
 115 ! Check if the residues involved in bridges are in the specified list of
 116 ! bridging residues.
 117         do i=1,nss
 118           do j=1,i-1
 119             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
 120             .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 121               write (iout,'(a,i3,a)') 'Disulfide pair',i,&
 122             ' contains residues present in other pairs.'
 123               write (*,'(a,i3,a)') 'Disulfide pair',i,&
 124             ' contains residues present in other pairs.'
 125 #ifdef MPI
 126               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 127               stop
 128 #endif
 129             endif
 130           enddo
 131           do j=1,ns
 132             if (ihpb(i).eq.iss(j)) goto 10
 133           enddo
 134           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 135    10     continue
 136           do j=1,ns
 137             if (jhpb(i).eq.iss(j)) goto 20
 138           enddo
 139           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 140    20     continue
 141           dhpb(i)=d0cm
 142           forcon(i)=akcm
 143         enddo
 144         do i=1,nss
 145           ihpb(i)=ihpb(i)+nres
 146           jhpb(i)=jhpb(i)+nres
 147         enddo
 148       endif
 149       endif
 150       write(iout,*) "end read_bridge"
 151       return
 152       end subroutine read_bridge
 153 !-----------------------------------------------------------------------------
 154       subroutine read_x(kanal,*)
 155
 156      use geometry_data
 157      use energy_data
 158      use geometry, only:int_from_cart1
 159 !      implicit real*8 (a-h,o-z)
 160 !      include 'DIMENSIONS'
 161 !      include 'COMMON.GEO'
 162 !      include 'COMMON.VAR'
 163 !      include 'COMMON.CHAIN'
 164 !      include 'COMMON.IOUNITS'
 165 !      include 'COMMON.CONTROL'
 166 !      include 'COMMON.LOCAL'
 167 !      include 'COMMON.INTERACT'
 168 ! Read coordinates from input
 169 !
 170 !el local variables
 171       integer :: l,k,j,i,kanal
 172
 173       read(kanal,'(8f10.5)',end=10,err=10) &
 174         ((c(l,k),l=1,3),k=1,nres),&
 175         ((c(l,k+nres),l=1,3),k=nnt,nct)
 176       do j=1,3
 177         c(j,nres+1)=c(j,1)
 178         c(j,2*nres)=c(j,nres)
 179       enddo
 180       call int_from_cart1(.false.)
 181       do i=1,nres-1
 182         do j=1,3
 183           dc(j,i)=c(j,i+1)-c(j,i)
 184           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 185         enddo
 186       enddo
 187       do i=nnt,nct
 188         if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
 189           do j=1,3
 190             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 191             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 192           enddo
 193         endif
 194       enddo
 195
 196       return
 197    10 return 1
 198       end subroutine read_x
 199 !-----------------------------------------------------------------------------
 200       subroutine read_threadbase
 201
 202       use geometry_data
 203       use energy_data
 204       use compare_data
 205 !      implicit real*8 (a-h,o-z)
 206 !      include 'DIMENSIONS'
 207 !      include 'COMMON.IOUNITS'
 208 !      include 'COMMON.GEO'
 209 !      include 'COMMON.VAR'
 210 !      include 'COMMON.INTERACT'
 211 !      include 'COMMON.LOCAL'
 212 !      include 'COMMON.NAMES'
 213 !      include 'COMMON.CHAIN'
 214 !      include 'COMMON.FFIELD'
 215 !      include 'COMMON.SBRIDGE'
 216 !      include 'COMMON.HEADER'
 217 !      include 'COMMON.CONTROL'
 218 !      include 'COMMON.DBASE'
 219 !      include 'COMMON.THREAD'
 220 !      include 'COMMON.TIME1'
 221
 222 !el local variables
 223       integer :: k,j,i
 224
 225 ! Read pattern database for threading.
 226       read (icbase,*) nseq
 227       allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
 228       allocate(nres_base(3,nseq)) !(3,maxseq)
 229       allocate(str_nam(nseq)) !(maxseq)
 230       do i=1,nseq
 231         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
 232          nres_base(2,i),nres_base(3,i)
 233         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
 234          nres_base(1,i))
 235 !       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
 236 !    &   nres_base(2,i),nres_base(3,i)
 237 !       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
 238 !    &   nres_base(1,i))
 239       enddo
 240       close (icbase)
 241       if (weidis.eq.0.0D0) weidis=0.1D0
 242       do i=nnt,nct
 243         do j=i+2,nct
 244           nhpb=nhpb+1
 245           ihpb(nhpb)=i
 246           jhpb(nhpb)=j
 247           forcon(nhpb)=weidis
 248         enddo
 249       enddo
 250       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
 251       write (iout,'(a,i5)') 'nexcl: ',nexcl
 252       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
 253       return
 254       end subroutine read_threadbase
 255 !-----------------------------------------------------------------------------
 256 #ifdef WHAM_RUN
 257 !el      subroutine read_angles(kanal,iscor,energ,iprot,*)
 258       subroutine read_angles(kanal,*)
 259
 260       use geometry_data
 261       use energy_data
 262 !      implicit real*8 (a-h,o-z)
 263 !      include 'DIMENSIONS'
 264 !      include 'DIMENSIONS.ZSCOPT'
 265 !      include 'COMMON.INTERACT'
 266 !      include 'COMMON.SBRIDGE'
 267 !      include 'COMMON.GEO'
 268 !      include 'COMMON.VAR'
 269 !      include 'COMMON.CHAIN'
 270 !      include 'COMMON.IOUNITS'
 271       character(len=80) :: lineh
 272       integer :: iscor,iprot,ic
 273       real(kind=8) :: energ
 274 #else
 275       subroutine read_angles(kanal,*)
 276
 277       use geometry_data
 278    !  use energy
 279    !  use control
 280 !      implicit real*8 (a-h,o-z)
 281 !      include 'DIMENSIONS'
 282 !      include 'COMMON.GEO'
 283 !      include 'COMMON.VAR'
 284 !      include 'COMMON.CHAIN'
 285 !      include 'COMMON.IOUNITS'
 286 !      include 'COMMON.CONTROL'
 287 #endif
 288 ! Read angles from input
 289 !
 290 !el local variables
 291       integer :: i,kanal
 292 #ifdef WHAM_RUN
 293       read(kanal,'(a80)',end=10,err=10) lineh
 294       read(lineh(:5),*,err=8) ic
 295       read(lineh(6:),*,err=8) energ
 296       goto 9
 297     8 ic=1
 298       print *,'error, assuming e=1d10',lineh
 299       energ=1d10
 300       nss=0
 301     9 continue
 302       read(lineh(18:),*,end=10,err=10) nss
 303       IF (NSS.LT.9) THEN
 304         read (lineh(20:),*,end=10,err=10) &
 305         (IHPB(I),JHPB(I),I=1,NSS),iscor
 306       ELSE
 307         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 308         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
 309           I=9,NSS),iscor
 310       ENDIF
 311 !      print *,"energy",energ," iscor",iscor
 312 #endif
 313        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 314        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 315        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 316        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 317
 318        do i=1,nres
 319 ! 9/7/01 avoid 180 deg valence angle
 320         if (theta(i).gt.179.99d0) theta(i)=179.99d0
 321 !
 322         theta(i)=deg2rad*theta(i)
 323         phi(i)=deg2rad*phi(i)
 324         alph(i)=deg2rad*alph(i)
 325         omeg(i)=deg2rad*omeg(i)
 326        enddo
 327       return
 328    10 return 1
 329       end subroutine read_angles
 330 !-----------------------------------------------------------------------------
 331       subroutine reada(rekord,lancuch,wartosc,default)
 332
 333 !      implicit none
 334       character*(*) :: rekord,lancuch
 335       real(kind=8) :: wartosc,default
 336       integer :: iread !,ilen
 337 !el      external ilen
 338       iread=index(rekord,lancuch)
 339       if (iread.eq.0) then
 340         wartosc=default
 341         return
 342       endif
 343       iread=iread+ilen(lancuch)+1
 344       read (rekord(iread:),*,err=10,end=10) wartosc
 345       return
 346   10  wartosc=default
 347       return
 348       end subroutine reada
 349 !-----------------------------------------------------------------------------
 350       subroutine readi(rekord,lancuch,wartosc,default)
 351
 352 !      implicit none
 353       character*(*) :: rekord,lancuch
 354       integer :: wartosc,default
 355       integer :: iread !,ilen
 356 !el      external ilen
 357       iread=index(rekord,lancuch)
 358       if (iread.eq.0) then
 359         wartosc=default
 360         return
 361       endif
 362       iread=iread+ilen(lancuch)+1
 363       read (rekord(iread:),*,err=10,end=10) wartosc
 364       return
 365   10  wartosc=default
 366       return
 367       end subroutine readi
 368 !-----------------------------------------------------------------------------
 369       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 370
 371 !      implicit none
 372       integer :: dim,i
 373       integer :: tablica(dim),default
 374       character*(*) :: rekord,lancuch
 375       character(len=80) :: aux
 376       integer :: iread !,ilen
 377 !el      external ilen
 378       do i=1,dim
 379         tablica(i)=default
 380       enddo
 381       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 382       if (iread.eq.0) return
 383       iread=iread+ilen(lancuch)+1
 384       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 385    10 return
 386       end subroutine multreadi
 387 !-----------------------------------------------------------------------------
 388       subroutine multreada(rekord,lancuch,tablica,dim,default)
 389
 390 !      implicit none
 391       integer :: dim,i
 392       real(kind=8) :: tablica(dim),default
 393       character*(*) :: rekord,lancuch
 394       character(len=80) :: aux
 395       integer :: iread !,ilen
 396 !el      external ilen
 397       do i=1,dim
 398         tablica(i)=default
 399       enddo
 400       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 401       if (iread.eq.0) return
 402       iread=iread+ilen(lancuch)+1
 403       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 404    10 return
 405       end subroutine multreada
 406 !-----------------------------------------------------------------------------
 407       subroutine card_concat(card,to_upper)
 408
 409 ! dla UNRESA to_upper jest zawsze .true.
 410 !      implicit real*8 (a-h,o-z)
 411 !      include 'DIMENSIONS'
 412 !      include 'COMMON.IOUNITS'
 413       character(*) :: card
 414       character(len=80) :: karta        !,ucase
 415       logical :: to_upper
 416 !el      external ilen
 417       read (inp,'(a)') karta
 418       if (to_upper) karta=ucase(karta)
 419       card=' '
 420       do while (karta(80:80).eq.'&')
 421         card=card(:ilen(card)+1)//karta(:79)
 422         read (inp,'(a)') karta
 423         if (to_upper) karta=ucase(karta)
 424       enddo
 425       card=card(:ilen(card)+1)//karta
 426       return
 427       end subroutine card_concat
 428 !----------------------------------------------------------------------------
 429       subroutine read_afminp
 430       use geometry_data
 431       use energy_data
 432       use control_data, only:out1file
 433       use MPI_data
 434       character*320 afmcard
 435       integer i
 436       print *, "wchodze"
 437       call card_concat(afmcard,.true.)
 438       call readi(afmcard,"BEG",afmbeg,0)
 439       call readi(afmcard,"END",afmend,0)
 440       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
 441       call reada(afmcard,"VEL",velAFMconst,0.0d0)
 442       print *,'FORCE=' ,forceAFMconst
 443 !------ NOW PROPERTIES FOR AFM
 444        distafminit=0.0d0
 445        do i=1,3
 446         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
 447        enddo
 448         distafminit=dsqrt(distafminit)
 449         print *,'initdist',distafminit
 450       return
 451       end subroutine read_afminp
 452 !-----------------------------------------------------------------------------
 453       subroutine read_dist_constr
 454       use MPI_data
 455    !  use control
 456       use geometry, only: dist
 457       use geometry_data
 458       use control_data
 459       use energy_data
 460 !      implicit real*8 (a-h,o-z)
 461 !      include 'DIMENSIONS'
 462 #ifdef MPI
 463       include 'mpif.h'
 464 #endif
 465 !      include 'COMMON.SETUP'
 466 !      include 'COMMON.CONTROL'
 467 !      include 'COMMON.CHAIN'
 468 !      include 'COMMON.IOUNITS'
 469 !      include 'COMMON.SBRIDGE'
 470       integer,dimension(2,100) :: ifrag_,ipair_
 471       real(kind=8),dimension(100) :: wfrag_,wpair_
 472       character(len=640) :: controlcard
 473
 474 !el local variables
 475       integer :: i,k,j,ddjk,ii,jj,itemp
 476       integer :: nfrag_,npair_,ndist_
 477       real(kind=8) :: dist_cut
 478
 479 !      write (iout,*) "Calling read_dist_constr"
 480 !      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 481 !      call flush(iout)
 482       if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
 483       if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
 484       if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
 485       if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
 486       if(.not.allocated(forcon)) allocate(forcon(maxdim))
 487       if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
 488       if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
 489       call gen_dist_constr2
 490       go to 1712
 491       endif
 492       call card_concat(controlcard,.true.)
 493       call readi(controlcard,"NFRAG",nfrag_,0)
 494       call readi(controlcard,"NPAIR",npair_,0)
 495       call readi(controlcard,"NDIST",ndist_,0)
 496       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 497       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 498       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 499       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 500       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 501 !      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 502 !      write (iout,*) "IFRAG"
 503 !      do i=1,nfrag_
 504 !        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 505 !      enddo
 506 !      write (iout,*) "IPAIR"
 507 !      do i=1,npair_
 508 !        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 509 !      enddo
 510 !      if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
 511 !      if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
 512 !      if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
 513 !      if(.not.allocated(forcon)) allocate(forcon(maxdim))
 514
 515       call flush(iout)
 516       do i=1,nfrag_
 517         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 518         if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
 519           ifrag_(2,i)=nstart_sup+nsup-1
 520 !        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 521         call flush(iout)
 522         if (wfrag_(i).gt.0.0d0) then
 523         do j=ifrag_(1,i),ifrag_(2,i)-1
 524           do k=j+1,ifrag_(2,i)
 525 !            write (iout,*) "j",j," k",k
 526             ddjk=dist(j,k)
 527             if (constr_dist.eq.1) then
 528               nhpb=nhpb+1
 529               ihpb(nhpb)=j
 530               jhpb(nhpb)=k
 531               dhpb(nhpb)=ddjk
 532               forcon(nhpb)=wfrag_(i)
 533             else if (constr_dist.eq.2) then
 534               if (ddjk.le.dist_cut) then
 535                 nhpb=nhpb+1
 536                 ihpb(nhpb)=j
 537                 jhpb(nhpb)=k
 538                 dhpb(nhpb)=ddjk
 539                 forcon(nhpb)=wfrag_(i)
 540               endif
 541             else
 542               nhpb=nhpb+1
 543               ihpb(nhpb)=j
 544               jhpb(nhpb)=k
 545               dhpb(nhpb)=ddjk
 546               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 547             endif
 548 #ifdef MPI
 549             if (.not.out1file .or. me.eq.king) &
 550             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
 551              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 552 #else
 553             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
 554              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 555 #endif
 556           enddo
 557         enddo
 558         endif
 559       enddo
 560       do i=1,npair_
 561         if (wpair_(i).gt.0.0d0) then
 562         ii = ipair_(1,i)
 563         jj = ipair_(2,i)
 564         if (ii.gt.jj) then
 565           itemp=ii
 566           ii=jj
 567           jj=itemp
 568         endif
 569         do j=ifrag_(1,ii),ifrag_(2,ii)
 570           do k=ifrag_(1,jj),ifrag_(2,jj)
 571             nhpb=nhpb+1
 572             ihpb(nhpb)=j
 573             jhpb(nhpb)=k
 574             forcon(nhpb)=wpair_(i)
 575             dhpb(nhpb)=dist(j,k)
 576 #ifdef MPI
 577             if (.not.out1file .or. me.eq.king) &
 578             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 579              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 580 #else
 581             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 582              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 583 #endif
 584           enddo
 585         enddo
 586         endif
 587       enddo
 588       do i=1,ndist_
 589 !        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 590 !        if (forcon(nhpb+1).gt.0.0d0) then
 591 !          nhpb=nhpb+1
 592 !          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 593         if (constr_dist.eq.11) then
 594         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
 595           ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 596         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 597         else
 598 !C        print *,"in else"
 599         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
 600           ibecarb(i),forcon(nhpb+1)
 601         endif
 602         if (forcon(nhpb+1).gt.0.0d0) then
 603           nhpb=nhpb+1
 604           if (ibecarb(i).gt.0) then
 605             ihpb(i)=ihpb(i)+nres
 606             jhpb(i)=jhpb(i)+nres
 607           endif
 608           if (dhpb(nhpb).eq.0.0d0) &
 609             dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 610         endif
 611
 612 #ifdef MPI
 613           if (.not.out1file .or. me.eq.king) &
 614           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 615            nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 616 #else
 617           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 618            nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 619 #endif
 620       enddo
 621  1712 continue
 622       call flush(iout)
 623       return
 624       end subroutine read_dist_constr
 625 !-----------------------------------------------------------------------------
 626       subroutine gen_dist_constr2
 627       use MPI_data
 628    !  use control
 629       use geometry, only: dist
 630       use geometry_data
 631       use control_data
 632       use energy_data
 633       integer :: i,j
 634       real(kind=8) :: distance
 635       if (constr_dist.eq.11) then
 636              do i=nstart_sup,nstart_sup+nsup-1
 637               do j=i+2,nstart_sup+nsup-1
 638                  distance=dist(i,j)
 639                  if (distance.le.15.0) then
 640                  nhpb=nhpb+1
 641                  ihpb(nhpb)=i+nstart_seq-nstart_sup
 642                  jhpb(nhpb)=j+nstart_seq-nstart_sup
 643                  forcon(nhpb)=sqrt(0.04*distance)
 644                  fordepth(nhpb)=sqrt(40.0/distance)
 645                  dhpb(nhpb)=distance-0.1d0
 646                  dhpb1(nhpb)=distance+0.1d0
 647
 648 #ifdef MPI
 649           if (.not.out1file .or. me.eq.king) &
 650           write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
 651           nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 652 #else
 653           write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
 654           nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 655 #endif
 656             endif
 657              enddo
 658            enddo
 659       else
 660       do i=nstart_sup,nstart_sup+nsup-1
 661         do j=i+2,nstart_sup+nsup-1
 662           nhpb=nhpb+1
 663           ihpb(nhpb)=i+nstart_seq-nstart_sup
 664           jhpb(nhpb)=j+nstart_seq-nstart_sup
 665           forcon(nhpb)=weidis
 666           dhpb(nhpb)=dist(i,j)
 667         enddo
 668       enddo
 669       endif
 670       return
 671       end subroutine gen_dist_constr2
 672
 673 !-----------------------------------------------------------------------------
 674 #ifdef WINIFL
 675       subroutine flush(iu)
 676       return
 677       end subroutine flush
 678 #endif
 679 #ifdef AIX
 680       subroutine flush(iu)
 681       call flush_(iu)
 682       return
 683       end subroutine flush
 684 #endif
 685 !-----------------------------------------------------------------------------
 686       subroutine copy_to_tmp(source)
 687
 688 !      include "DIMENSIONS"
 689 !      include "COMMON.IOUNITS"
 690       character*(*) :: source
 691       character(len=256) :: tmpfile
 692 !      integer ilen
 693 !el      external ilen
 694       logical :: ex
 695       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
 696       inquire(file=tmpfile,exist=ex)
 697       if (ex) then
 698         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
 699          " to temporary directory..."
 700         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
 701         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
 702       endif
 703       return
 704       end subroutine copy_to_tmp
 705 !-----------------------------------------------------------------------------
 706       subroutine move_from_tmp(source)
 707
 708 !      include "DIMENSIONS"
 709 !      include "COMMON.IOUNITS"
 710       character*(*) :: source
 711 !      integer ilen
 712 !el      external ilen
 713       write (*,*) "Moving ",source(:ilen(source)),&
 714        " from temporary directory to working directory"
 715       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
 716       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
 717       return
 718       end subroutine move_from_tmp
 719 !-----------------------------------------------------------------------------
 720 ! misc.f
 721 !-----------------------------------------------------------------------------
 722 ! $Date: 1994/10/12 17:24:21 $
 723 ! $Revision: 2.5 $
 724
 725       logical function find_arg(ipos,line,errflag)
 726
 727       integer, parameter :: maxlen = 80
 728       character(len=80) :: line
 729       character(len=1) :: empty=' ',equal='='
 730       logical :: errflag
 731       integer :: ipos
 732 ! This function returns .TRUE., if an argument follows keyword keywd; if so
 733 ! IPOS will point to the first non-blank character of the argument. Returns
 734 ! .FALSE., if no argument follows the keyword; in this case IPOS points
 735 ! to the first non-blank character of the next keyword.
 736
 737       do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
 738         ipos=ipos+1
 739       enddo
 740       errflag=.false.
 741       if (line(ipos:ipos).eq.equal) then
 742          find_arg=.true.
 743          ipos=ipos+1
 744          do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
 745            ipos=ipos+1
 746          enddo
 747          if (ipos.gt.maxlen) errflag=.true.
 748       else
 749          find_arg=.false.
 750       endif
 751
 752       return
 753       end function find_arg
 754 !-----------------------------------------------------------------------------
 755       logical function find_group(iunit,jout,key1)
 756
 757       character*(*) :: key1
 758       character(len=80) :: karta        !,ucase
 759       integer :: iunit,jout
 760       integer :: ll !,ilen
 761 !EL      external ilen
 762 !EL      logical lcom
 763       rewind (iunit)
 764       karta=' '
 765       ll=ilen(key1)
 766       do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
 767         read (iunit,'(a)',end=10) karta
 768       enddo
 769       write (jout,'(2a)') '> ',karta(1:78)
 770       find_group=.true.
 771       return
 772    10 find_group=.false.
 773       return
 774       end function find_group
 775 !-----------------------------------------------------------------------------
 776       logical function iblnk(charc)
 777       character(len=1) :: charc
 778       integer :: n
 779       n = ichar(charc)
 780       iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
 781       return
 782       end function iblnk
 783 !-----------------------------------------------------------------------------
 784       integer function ilen(string)
 785       character*(*) ::  string
 786 !EL      logical :: iblnk
 787
 788       ilen = len(string)
 789 1     if ( ilen .gt. 0 ) then
 790          if ( iblnk( string(ilen:ilen) ) ) then
 791             ilen = ilen - 1
 792             goto 1
 793          endif
 794       endif
 795       return
 796       end function ilen
 797 !-----------------------------------------------------------------------------
 798       integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
 799       integer :: nkey,i,ikey,narg
 800       character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
 801 !      character(len=16) :: ucase
 802
 803       do i=1,narg
 804         if (ucase(keywd).eq.keywdset(i,ikey)) then
 805 ! Match found
 806           in_keywd_set=i
 807           return
 808         endif
 809       enddo
 810 ! No match to the allowed set of keywords if this point is reached.
 811       in_keywd_set=0
 812       return
 813       end function in_keywd_set
 814 !-----------------------------------------------------------------------------
 815       character function lcase(string)
 816       integer :: i, k, idiff
 817       character*(*) :: string
 818       character(len=1) :: c
 819       character(len=40) :: chtmp
 820 !
 821       i = len(lcase)
 822       k = len(string)
 823       if (i .lt. k) then
 824          k = i
 825          if (string(k+1:) .ne. ' ') then
 826             chtmp = string
 827          endif
 828       endif
 829       idiff = ichar('a') - ichar('A')
 830       lcase = string
 831       do 99 i = 1, k
 832          c = string(i:i)
 833          if (lge(c,'A') .and. lle(c,'Z')) then
 834             lcase(i:i) = char(ichar(c) + idiff)
 835          endif
 836    99 continue
 837       return
 838       end function lcase
 839 !-----------------------------------------------------------------------------
 840       logical function lcom(ipos,karta)
 841       character(len=80) :: karta
 842       character :: koment(2) = (/'!','#'/)
 843       integer :: ipos,i
 844
 845       lcom=.false.
 846       do i=1,2
 847         if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
 848       enddo
 849       return
 850       end function lcom
 851 !-----------------------------------------------------------------------------
 852       logical function lower_case(ch)
 853       character*(*) :: ch
 854       lower_case=(ch.ge.'a' .and. ch.le.'z')
 855       return
 856       end function lower_case
 857 !-----------------------------------------------------------------------------
 858       subroutine mykey(line,keywd,ipos,blankline,errflag)
 859
 860 ! This subroutine seeks a non-empty substring keywd in the string LINE.
 861 ! The substring begins with the first character different from blank and
 862 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
 863 ! at the character left to the first blank or "=". When the subroutine is
 864 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
 865 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
 866 ! only separators or the maximum length of the data line (80) has been reached.
 867 ! The logical variable ERRFLAG is set at .TRUE. if the string
 868 ! consists only from a "=".
 869       integer, parameter :: maxlen=80
 870       character(len=1) :: empty=' ',equal='=',comma=','
 871       character*(*) :: keywd
 872       character(len=80) :: line
 873       logical :: blankline,errflag      !EL,lcom
 874       integer :: ipos,istart,iend
 875       errflag=.false.
 876       do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
 877         ipos=ipos+1
 878       enddo
 879       if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
 880 ! At this point the rest of the input line turned out to contain only blanks
 881 ! or to be commented out.
 882         blankline=.true.
 883         return
 884       endif
 885       blankline=.false.
 886       istart=ipos
 887 ! Checks whether the current char is a separator.
 888       do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
 889        .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
 890         ipos=ipos+1
 891       enddo
 892       iend=ipos-1
 893 ! Error flag set to .true., if the length of the keyword was found less than 1.
 894       if (iend.lt.istart) then
 895         errflag=.true.
 896         return
 897       endif
 898       keywd=line(istart:iend)
 899       return
 900       end subroutine  mykey
 901 !-----------------------------------------------------------------------------
 902       subroutine numstr(inum,numm)
 903       character(len=10) :: huj='0123456789'
 904       character*(*) :: numm
 905       integer :: inumm,inum,inum1,inum2
 906       inumm=inum
 907       inum1=inumm/10
 908       inum2=inumm-10*inum1
 909       inumm=inum1
 910       numm(3:3)=huj(inum2+1:inum2+1)
 911       inum1=inumm/10
 912       inum2=inumm-10*inum1
 913       inumm=inum1
 914       numm(2:2)=huj(inum2+1:inum2+1)
 915       inum1=inumm/10
 916       inum2=inumm-10*inum1
 917       inumm=inum1
 918       numm(1:1)=huj(inum2+1:inum2+1)
 919       return
 920       end subroutine numstr
 921 !-----------------------------------------------------------------------------
 922       function ucase(string)
 923       integer :: i, k, idiff
 924       character(*) :: string
 925       character(len=len(string)) :: ucase
 926       character(len=1) :: c
 927       character(len=40) :: chtmp
 928 !
 929       i = len(ucase)
 930       k = len(string)
 931       if (i .lt. k) then
 932          k = i
 933          if (string(k+1:) .ne. ' ') then
 934             chtmp = string
 935          endif
 936       endif
 937       idiff = ichar('a') - ichar('A')
 938       ucase = string
 939       do 99 i = 1, k
 940          c = string(i:i)
 941          if (lge(c,'a') .and. lle(c,'z')) then
 942             ucase(i:i) = char(ichar(c) - idiff)
 943          endif
 944    99 continue
 945       return
 946       end function ucase
 947 !-----------------------------------------------------------------------------
 948 ! geomout.F
 949 !-----------------------------------------------------------------------------
 950       subroutine pdbout(etot,tytul,iunit)
 951
 952       use geometry_data, only: c,nres
 953       use energy_data
 954    !  use control
 955       use compare_data
 956       use MD_data
 957 !      implicit real*8 (a-h,o-z)
 958 !      include 'DIMENSIONS'
 959 !      include 'COMMON.CHAIN'
 960 !      include 'COMMON.INTERACT'
 961 !      include 'COMMON.NAMES'
 962 !      include 'COMMON.IOUNITS'
 963 !      include 'COMMON.HEADER'
 964 !      include 'COMMON.SBRIDGE'
 965 !      include 'COMMON.DISTFIT'
 966 !      include 'COMMON.MD'
 967 !el      character(len=50) :: tytul
 968       character*(*) :: tytul
 969       character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
 970       integer,dimension(nres) :: ica    !(maxres)
 971        integer iti1
 972 !el  local variables
 973       integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
 974       real(kind=8) :: etot
 975       integer :: nres2
 976       nres2=2*nres
 977
 978       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
 979
 980       write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
 981 !model      write (iunit,'(a5,i6)') 'MODEL',1
 982       if (nhfrag.gt.0) then
 983        do j=1,nhfrag
 984         iti=itype(hfrag(1,j),1)
 985         itj=itype(hfrag(2,j),1)
 986         if (j.lt.10) then
 987            write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
 988                  'HELIX',j,'H',j,&
 989                  restyp(iti,1),hfrag(1,j)-1,&
 990                  restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
 991         else
 992              write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
 993                  'HELIX',j,'H',j,&
 994                  restyp(iti,1),hfrag(1,j)-1,&
 995                  restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
 996         endif
 997        enddo
 998       endif
 999
1000       if (nbfrag.gt.0) then
1001
1002        do j=1,nbfrag
1003
1004         iti=itype(bfrag(1,j),1)
1005         itj=itype(bfrag(2,j)-1,1)
1006
1007         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1008                  'SHEET',1,'B',j,2,&
1009                  restyp(iti,1),bfrag(1,j)-1,&
1010                  restyp(itj,1),bfrag(2,j)-2,0
1011
1012         if (bfrag(3,j).gt.bfrag(4,j)) then
1013
1014          itk=itype(bfrag(3,j),1)
1015          itl=itype(bfrag(4,j)+1,1)
1016
1017          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)') &
1018                  'SHEET',2,'B',j,2,&
1019                  restyp(itl,1),bfrag(4,j),&
1020                  restyp(itk,1),bfrag(3,j)-1,-1,&
1021                  "N",restyp(itk,1),bfrag(3,j)-1,&
1022                  "O",restyp(iti,1),bfrag(1,j)-1
1023
1024         else
1025
1026          itk=itype(bfrag(3,j),1)
1027          itl=itype(bfrag(4,j)-1,1)
1028
1029
1030         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)') &
1031                  'SHEET',2,'B',j,2,&
1032                  restyp(itk,1),bfrag(3,j)-1,&
1033                  restyp(itl,1),bfrag(4,j)-2,1,&
1034                  "N",restyp(itk,1),bfrag(3,j)-1,&
1035                  "O",restyp(iti,1),bfrag(1,j)-1
1036
1037
1038
1039         endif
1040
1041        enddo
1042       endif
1043
1044       if (nss.gt.0) then
1045         do i=1,nss
1046          if (dyn_ss) then
1047           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1048                'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1049                           'CYS',jdssb(i)-nnt+1
1050          else
1051           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1052                'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1053                           'CYS',jhpb(i)-nnt+1-nres
1054          endif
1055         enddo
1056       endif
1057
1058       iatom=0
1059       ichain=1
1060       ires=0
1061       do i=nnt,nct
1062         iti=itype(i,molnum(i))
1063         print *,i,molnum(i)
1064         if (molnum(i+1).eq.0) then
1065         iti1=ntyp1_molec(molnum(i))
1066         else
1067         iti1=itype(i+1,molnum(i+1))
1068         endif
1069         if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle
1070         if (i.lt.nnt) then
1071         if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle
1072         endif
1073         if (iti.eq.ntyp1) then
1074           ichain=ichain+1
1075 !          ires=0
1076           write (iunit,'(a)') 'TER'
1077         else
1078         ires=ires+1
1079         iatom=iatom+1
1080         ica(i)=iatom
1081         if (molnum(i).ne.5) then
1082         write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1083            ires,(c(j,i),j=1,3),vtot(i)
1084         else
1085         write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1086            ires,(c(j,i),j=1,3),vtot(i)
1087         endif
1088         if ((iti.ne.10).and.(molnum(i).ne.5)) then
1089           iatom=iatom+1
1090           write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1091             ires,(c(j,nres+i),j=1,3),&
1092             vtot(i+nres)
1093         endif
1094         endif
1095       enddo
1096       write (iunit,'(a)') 'TER'
1097       do i=nnt,nct-1
1098         if (itype(i,1).eq.ntyp1) cycle
1099         if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1100           write (iunit,30) ica(i),ica(i+1)
1101         else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1102           write (iunit,30) ica(i),ica(i+1),ica(i)+1
1103         else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1104           write (iunit,30) ica(i),ica(i)+1
1105         endif
1106       enddo
1107       if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1108         write (iunit,30) ica(nct),ica(nct)+1
1109       endif
1110       do i=1,nss
1111        if (dyn_ss) then
1112         write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1113        else
1114         write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1115        endif
1116       enddo
1117       write (iunit,'(a6)') 'ENDMDL'
1118   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1119
1120   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1121   30  FORMAT ('CONECT',8I5)
1122   60  FORMAT ('HETATM',I5,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1123
1124       return
1125       end subroutine pdbout
1126 !-----------------------------------------------------------------------------
1127       subroutine MOL2out(etot,tytul)
1128 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1129 ! format.
1130       use geometry_data, only: c
1131       use energy_data
1132    !  use control
1133 !      implicit real*8 (a-h,o-z)
1134 !      include 'DIMENSIONS'
1135 !      include 'COMMON.CHAIN'
1136 !      include 'COMMON.INTERACT'
1137 !      include 'COMMON.NAMES'
1138 !      include 'COMMON.IOUNITS'
1139 !      include 'COMMON.HEADER'
1140 !      include 'COMMON.SBRIDGE'
1141       character(len=32) :: tytul,fd
1142       character(len=3) :: zahl
1143       character(len=6) :: res_num,pom   !,ucase
1144
1145 !el  local variables
1146       integer :: i,j
1147       real(kind=8) :: etot
1148
1149 #ifdef AIX
1150       call fdate_(fd)
1151 #elif (defined CRAY)
1152       call date(fd)
1153 #else
1154       call fdate(fd)
1155 #endif
1156       write (imol2,'(a)') '#'
1157       write (imol2,'(a)') &
1158        '#         Creating user name:           unres'
1159       write (imol2,'(2a)') '#         Creation time:                ',&
1160        fd
1161       write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1162       write (imol2,'(a)') tytul
1163       write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1164       write (imol2,'(a)') 'SMALL'
1165       write (imol2,'(a)') 'USER_CHARGES'
1166       write (imol2,'(a)') '\@<TRIPOS>ATOM'
1167       do i=nnt,nct
1168         write (zahl,'(i3)') i
1169         pom=ucase(restyp(itype(i,1),1))
1170         res_num = pom(:3)//zahl(2:)
1171         write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1172       enddo
1173       write (imol2,'(a)') '\@<TRIPOS>BOND'
1174       do i=nnt,nct-1
1175         write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1176       enddo
1177       do i=1,nss
1178         write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1179       enddo
1180       write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1181       do i=nnt,nct
1182         write (zahl,'(i3)') i
1183         pom = ucase(restyp(itype(i,1),1))
1184         res_num = pom(:3)//zahl(2:)
1185         write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1186       enddo
1187   10  FORMAT (I7,' CA      ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1188   30  FORMAT (I7,1x,A,I14,' RESIDUE',I13,' ****  ****')
1189       return
1190       end subroutine MOL2out
1191 !-----------------------------------------------------------------------------
1192       subroutine intout
1193
1194       use geometry_data
1195       use energy_data, only: itype
1196    !  use control
1197 !      implicit real*8 (a-h,o-z)
1198 !      include 'DIMENSIONS'
1199 !      include 'COMMON.IOUNITS'
1200 !      include 'COMMON.CHAIN'
1201 !      include 'COMMON.VAR'
1202 !      include 'COMMON.LOCAL'
1203 !      include 'COMMON.INTERACT'
1204 !      include 'COMMON.NAMES'
1205 !      include 'COMMON.GEO'
1206 !      include 'COMMON.TORSION'
1207 !el  local variables
1208       integer :: i,iti
1209
1210       write (iout,'(/a)') 'Geometry of the virtual chain.'
1211       write (iout,'(7a)') '  Res  ','         d','     Theta',&
1212        '       Phi','       Dsc','     Alpha','      Omega'
1213       do i=1,nres
1214         iti=itype(i,1)
1215 !         print *,vbld(i),"vbld(i)",i
1216         write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
1217            rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1218            rad2deg*omeg(i)
1219       enddo
1220       return
1221       end subroutine intout
1222 !-----------------------------------------------------------------------------
1223 #ifdef CLUSTER
1224       subroutine briefout(it,ener,free)!,plik)
1225 #else
1226       subroutine briefout(it,ener)
1227 #endif
1228       use geometry_data
1229       use energy_data
1230    !  use control
1231 !      implicit real*8 (a-h,o-z)
1232 !      include 'DIMENSIONS'
1233 !      include 'COMMON.IOUNITS'
1234 !      include 'COMMON.CHAIN'
1235 !      include 'COMMON.VAR'
1236 !      include 'COMMON.LOCAL'
1237 !      include 'COMMON.INTERACT'
1238 !      include 'COMMON.NAMES'
1239 !      include 'COMMON.GEO'
1240 !      include 'COMMON.SBRIDGE'
1241 !     print '(a,i5)',intname,igeom
1242 !el  local variables
1243       integer :: i,it
1244       real(kind=8) :: ener,free
1245 !     character(len=80) :: plik
1246 !      integer :: iii
1247
1248 #if defined(AIX) || defined(PGI)
1249       open (igeom,file=intname,position='append')
1250 #else
1251       open (igeom,file=intname,access='append')
1252 #endif
1253 #ifdef WHAM_RUN
1254 !      iii=igeom
1255       igeom=iout
1256 #endif
1257       print *,nss
1258       IF (NSS.LE.9) THEN
1259 #ifdef CLUSTER
1260         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1261       ELSE
1262         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1263 #else
1264         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1265       ELSE
1266         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1267 #endif
1268         WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1269       ENDIF
1270 !     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1271       WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1272       WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1273 !     if (nvar.gt.nphi+ntheta) then
1274         write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1275         write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1276 !     endif
1277       close(igeom)
1278   180 format (I5,F12.3,I2,9(1X,2I3))
1279   190 format (3X,11(1X,2I3))
1280   200 format (8F10.4)
1281       return
1282       end subroutine briefout
1283 !-----------------------------------------------------------------------------
1284 #ifdef WINIFL
1285       subroutine fdate(fd)
1286       character(len=32) :: fd
1287       write(fd,'(32x)')
1288       return
1289       end subroutine fdate
1290 #endif
1291 !-----------------------------------------------------------------------------
1292 #ifdef WHAM_RUN
1293       real(kind=8) function gyrate(jcon)
1294 #else
1295       real(kind=8) function gyrate()
1296 #endif
1297
1298       use geometry_data, only: c
1299    !  use geometry
1300       use energy_data
1301 !      implicit real*8 (a-h,o-z)
1302 !      include 'DIMENSIONS'
1303 !      include 'COMMON.INTERACT'
1304 !      include 'COMMON.CHAIN'
1305       real(kind=8) :: rg
1306       real(kind=8),dimension(3) :: cen
1307 !el  local variables
1308       integer :: i,j,jcon
1309
1310       do j=1,3
1311        cen(j)=0.0d0
1312       enddo
1313
1314       do i=nnt,nct
1315           do j=1,3
1316             cen(j)=cen(j)+c(j,i)
1317           enddo
1318       enddo
1319       do j=1,3
1320             cen(j)=cen(j)/dble(nct-nnt+1)
1321       enddo
1322       rg = 0.0d0
1323       do i = nnt, nct
1324         do j=1,3
1325          rg = rg + (c(j,i)-cen(j))**2
1326         enddo
1327       end do
1328 #ifdef WHAM_RUN
1329       gyrate = dsqrt(rg/dble(nct-nnt+1))
1330 #else
1331       gyrate = sqrt(rg/dble(nct-nnt+1))
1332 #endif
1333       return
1334       end function gyrate
1335 #ifdef WHAM_RUN
1336 !-----------------------------------------------------------------------------
1337 ! readrtns.F WHAM
1338       subroutine reads(rekord,lancuch,wartosc,default)
1339 !      implicit none
1340       character*(*) :: rekord,lancuch,wartosc,default
1341       character(len=80) :: aux
1342       integer :: lenlan,lenrec,iread,ireade
1343 !el      external ilen
1344 !el      logical iblnk
1345 !el      external iblnk
1346       lenlan=ilen(lancuch)
1347       lenrec=ilen(rekord)
1348       iread=index(rekord,lancuch(:lenlan)//"=")
1349 !      print *,"rekord",rekord," lancuch",lancuch
1350 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1351       if (iread.eq.0) then
1352         wartosc=default
1353         return
1354       endif
1355       iread=iread+lenlan+1
1356 !      print *,"iread",iread
1357 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1358       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1359         iread=iread+1
1360 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1361       enddo
1362 !      print *,"iread",iread
1363       if (iread.gt.lenrec) then
1364          wartosc=default
1365         return
1366       endif
1367       ireade=iread+1
1368 !      print *,"ireade",ireade
1369       do while (ireade.lt.lenrec .and. &
1370          .not.iblnk(rekord(ireade:ireade)))
1371         ireade=ireade+1
1372       enddo
1373       wartosc=rekord(iread:ireade)
1374       return
1375       end subroutine reads
1376 #endif
1377 !-----------------------------------------------------------------------------
1378 ! permut.F
1379 !-----------------------------------------------------------------------------
1380       subroutine permut(isym)
1381
1382       use geometry_data, only: tabperm
1383 !      implicit real*8 (a-h,o-z)
1384 !      include 'DIMENSIONS'
1385 !      include 'COMMON.LOCAL'
1386 !      include 'COMMON.VAR'
1387 !      include 'COMMON.CHAIN'
1388 !      include 'COMMON.INTERACT'
1389 !      include 'COMMON.IOUNITS'
1390 !      include 'COMMON.GEO'
1391 !      include 'COMMON.NAMES'
1392 !      include 'COMMON.CONTROL'
1393
1394       integer :: n,isym
1395 !      logical nextp
1396 !el      external nextp
1397       integer,dimension(isym) :: a
1398 !      parameter(n=symetr)
1399 !el local variables
1400       integer :: kkk,i
1401
1402       n=isym
1403       if (n.eq.1) then
1404         tabperm(1,1)=1
1405         return
1406       endif
1407       kkk=0
1408       do i=1,n
1409       a(i)=i
1410       enddo
1411    10 print *,(a(i),i=1,n)
1412       kkk=kkk+1
1413       do i=1,n
1414       tabperm(kkk,i)=a(i)
1415 !      write (iout,*) "tututu", kkk
1416       enddo
1417       if(nextp(n,a)) go to 10
1418       return
1419       end subroutine permut
1420 !-----------------------------------------------------------------------------
1421       logical function nextp(n,a)
1422
1423       integer :: n,i,j,k,t
1424 !      logical :: nextp
1425       integer,dimension(n) :: a
1426       i=n-1
1427    10 if(a(i).lt.a(i+1)) go to 20
1428       i=i-1
1429       if(i.eq.0) go to 20
1430       go to 10
1431    20 j=i+1
1432       k=n
1433    30 t=a(j)
1434       a(j)=a(k)
1435       a(k)=t
1436       j=j+1
1437       k=k-1
1438       if(j.lt.k) go to 30
1439       j=i
1440       if(j.ne.0) go to 40
1441       nextp=.false.
1442       return
1443    40 j=j+1
1444       if(a(j).lt.a(i)) go to 40
1445       t=a(i)
1446       a(i)=a(j)
1447       a(j)=t
1448       nextp=.true.
1449       return
1450       end function nextp
1451 !-----------------------------------------------------------------------------
1452 !-----------------------------------------------------------------------------
1453       end module io_base