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