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