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