source/unres/io_base.f90

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