source/unres/io.F90

   1       module io
   2 !-----------------------------------------------------------------------
   3       use io_units
   4       use names
   5       use io_base
   6       use io_config
   7       implicit none
   8 !-----------------------------------------------------------------------------
   9 !
  10 !
  11 !-----------------------------------------------------------------------------
  12       contains
  13 !-----------------------------------------------------------------------------
  14 ! bank.F    io_csa
  15 !-----------------------------------------------------------------------------
  16       subroutine write_csa_pdb(var,ene,nft,ik,iw_pdb)
  17
  18       use csa_data
  19       use geometry_data, only:nres,nvar
  20       use geometry, only:var_to_geom,chainbuild
  21       use compare, only:secondary2
  22 !      implicit real*8 (a-h,o-z)
  23 !      include 'DIMENSIONS'
  24 !      include 'COMMON.CSA'
  25 !      include 'COMMON.BANK'
  26 !      include 'COMMON.VAR'
  27 !      include 'COMMON.IOUNITS'
  28 !      include 'COMMON.MINIM'
  29 !      include 'COMMON.SETUP'
  30 !      include 'COMMON.GEO'
  31 !      include 'COMMON.CHAIN'
  32 !      include 'COMMON.LOCAL'
  33 !      include 'COMMON.INTERACT'
  34 !      include 'COMMON.NAMES'
  35 !      include 'COMMON.SBRIDGE'
  36       integer :: lenpre,lenpot  !,ilen
  37 !el      external ilen
  38       real(kind=8),dimension(nvar) :: var       !(maxvar)       (maxvar=6*maxres)
  39       character(len=50) :: titelloc
  40       character(len=3) :: zahl
  41       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: ene
  42 !el local variables
  43       integer :: nft,ik,iw_pdb
  44
  45       nmin_csa=nmin_csa+1
  46       if(ene(1).lt.eglob_csa) then
  47         eglob_csa=ene(1)
  48         nglob_csa=nglob_csa+1
  49         call numstr(nglob_csa,zahl)
  50
  51         call var_to_geom(nvar,var)
  52         call chainbuild
  53         call secondary2(.false.)
  54
  55         lenpre=ilen(prefix)
  56         open(icsa_pdb,file=prefix(:lenpre)//'@'//zahl//'.pdb')
  57
  58         if (iw_pdb.eq.1) then
  59           write(titelloc,'(a2,i3,a3,i9,a3,i6)') &
  60           'GM',nglob_csa,' e ',nft,' m ',nmin_csa
  61         else
  62           write(titelloc,'(a2,i3,a3,i9,a3,i6,a5,f5.2,a5,f5.1)') &
  63          'GM',nglob_csa,' e ',nft,' m ',nmin_csa,' rms ',&
  64                rmsn(ik),' %NC ',pncn(ik)*100
  65         endif
  66         call pdbout(eglob_csa,titelloc,icsa_pdb)
  67         close(icsa_pdb)
  68       endif
  69
  70       return
  71       end subroutine write_csa_pdb
  72 !-----------------------------------------------------------------------------
  73 ! csa.f         io_csa
  74 !-----------------------------------------------------------------------------
  75       subroutine from_pdb(n,idum)
  76 ! This subroutine stores the UNRES int variables generated from
  77 ! subroutine readpdb into the 1st conformation of in dihang_in.
  78 ! Subsequent n-1 conformations of dihang_in have identical values
  79 ! of theta and phi as the 1st conformation but random values for
  80 ! alph and omeg.
  81 ! The array cref (also generated from subroutine readpdb) is stored
  82 ! to crefjlee to be used for rmsd calculation in CSA, if necessary.
  83
  84       use csa_data
  85       use geometry_data
  86       use random, only: ran1
  87 !      implicit real*8 (a-h,o-z)
  88 !      include 'DIMENSIONS'
  89 !      include 'COMMON.IOUNITS'
  90 !      include 'COMMON.CHAIN'
  91 !      include 'COMMON.VAR'
  92 !      include 'COMMON.BANK'
  93 !      include 'COMMON.GEO'
  94 !el local variables
  95       integer :: n,idum,m,i,j,k,kk,kkk
  96       real(kind=8) :: e
  97
  98       m=1
  99       do j=2,nres-1
 100         dihang_in(1,j,1,m)=theta(j+1)
 101         dihang_in(2,j,1,m)=phi(j+2)
 102         dihang_in(3,j,1,m)=alph(j)
 103         dihang_in(4,j,1,m)=omeg(j)
 104       enddo
 105       dihang_in(2,nres-1,1,k)=0.0d0
 106
 107       do m=2,n
 108        do k=2,nres-1
 109         dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
 110         dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
 111         if(dabs(dihang_in(3,k,1,1)).gt.1.d-6) then
 112          dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
 113          dihang_in(3,k,1,m)=dihang_in(3,k,1,m)*deg2rad
 114         endif
 115         if(dabs(dihang_in(4,k,1,1)).gt.1.d-6) then
 116          dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
 117          dihang_in(4,k,1,m)=dihang_in(4,k,1,m)*deg2rad
 118         endif
 119        enddo
 120       enddo
 121
 122 ! Store cref to crefjlee (they are in COMMON.CHAIN).
 123       do k=1,2*nres
 124        do kk=1,3
 125         kkk=1
 126         crefjlee(kk,k)=cref(kk,k,kkk)
 127        enddo
 128       enddo
 129
 130       open(icsa_native_int,file=csa_native_int,status="old")
 131       do m=1,n
 132          write(icsa_native_int,*) m,e
 133          write(icsa_native_int,200) &
 134           (dihang_in(1,k,1,m)*rad2deg,k=2,nres-1)
 135          write(icsa_native_int,200) &
 136           (dihang_in(2,k,1,m)*rad2deg,k=2,nres-2)
 137          write(icsa_native_int,200) &
 138           (dihang_in(3,k,1,m)*rad2deg,k=2,nres-1)
 139          write(icsa_native_int,200) &
 140           (dihang_in(4,k,1,m)*rad2deg,k=2,nres-1)
 141       enddo
 142
 143       do k=1,nres
 144        write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
 145       enddo
 146       close(icsa_native_int)
 147
 148   200 format (8f10.4)
 149
 150       return
 151       end subroutine from_pdb
 152 !-----------------------------------------------------------------------------
 153       subroutine from_int(n,mm,idum)
 154
 155       use csa_data
 156       use geometry_data
 157       use energy_data
 158       use geometry, only:chainbuild,gen_side
 159       use energy, only:etotal
 160       use compare
 161 !      implicit real*8 (a-h,o-z)
 162 !      include 'DIMENSIONS'
 163 !      include 'COMMON.IOUNITS'
 164 !      include 'COMMON.CHAIN'
 165 !      include 'COMMON.VAR'
 166 !      include 'COMMON.INTERACT'
 167 !      include 'COMMON.BANK'
 168 !      include 'COMMON.GEO'
 169 !      include 'COMMON.CONTACTS'
 170 !      integer ilen
 171 !el      external ilen
 172       logical :: fail
 173       real(kind=8),dimension(0:n_ene) :: energia
 174 !el local variables
 175       integer :: n,mm,idum,i,ii,j,m,k,kk,maxcount_fail,icount_fail,maxsi
 176       real(kind=8) :: co
 177
 178       open(icsa_native_int,file=csa_native_int,status="old")
 179       read (icsa_native_int,*)
 180       call read_angles(icsa_native_int,*10)
 181       goto 11
 182    10 write (iout,'(2a)') "CHUJ NASTAPIL - error in ",&
 183         csa_native_int(:ilen(csa_native_int))
 184    11 continue
 185       call intout
 186       do j=2,nres-1
 187         dihang_in(1,j,1,1)=theta(j+1)
 188         dihang_in(2,j,1,1)=phi(j+2)
 189         dihang_in(3,j,1,1)=alph(j)
 190         dihang_in(4,j,1,1)=omeg(j)
 191       enddo
 192       dihang_in(2,nres-1,1,1)=0.0d0
 193
 194 !         read(icsa_native_int,*) ind,e
 195 !         read(icsa_native_int,200) (dihang_in(1,k,1,1),k=2,nres-1)
 196 !         read(icsa_native_int,200) (dihang_in(2,k,1,1),k=2,nres-2)
 197 !         read(icsa_native_int,200) (dihang_in(3,k,1,1),k=2,nres-1)
 198 !         read(icsa_native_int,200) (dihang_in(4,k,1,1),k=2,nres-1)
 199 !         dihang_in(2,nres-1,1,1)=0.d0
 200
 201          maxsi=100
 202          maxcount_fail=100
 203
 204          do m=mm+2,n
 205 !          do k=2,nres-1
 206 !           dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
 207 !           dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
 208 !           if(abs(dihang_in(3,k,1,1)).gt.1.d-3) then
 209 !            dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
 210 !           endif
 211 !           if(abs(dihang_in(4,k,1,1)).gt.1.d-3) then
 212 !            dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
 213 !           endif
 214 !          enddo
 215 !           call intout
 216            fail=.true.
 217
 218            icount_fail=0
 219
 220            DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL)
 221
 222            do i=nnt,nct
 223              if (itype(i,1).ne.10) then
 224 !d             print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
 225                fail=.true.
 226                ii=0
 227                do while (fail .and. ii .le. maxsi)
 228                  call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail)
 229                  ii = ii+1
 230                enddo
 231              endif
 232            enddo
 233            call chainbuild
 234            call etotal(energia)
 235            fail = (energia(0).ge.1.0d20)
 236            icount_fail=icount_fail+1
 237
 238            ENDDO
 239
 240            if (icount_fail.gt.maxcount_fail) then
 241              write (iout,*) &
 242              'Failed to generate non-overlaping near-native conf.',&
 243              m
 244            endif
 245
 246            do j=2,nres-1
 247              dihang_in(1,j,1,m)=theta(j+1)
 248              dihang_in(2,j,1,m)=phi(j+2)
 249              dihang_in(3,j,1,m)=alph(j)
 250              dihang_in(4,j,1,m)=omeg(j)
 251            enddo
 252            dihang_in(2,nres-1,1,m)=0.0d0
 253          enddo
 254
 255 !      do m=1,n
 256 !        write(icsa_native_int,*) m,e
 257 !         write(icsa_native_int,200) (dihang_in(1,k,1,m),k=2,nres-1)
 258 !         write(icsa_native_int,200) (dihang_in(2,k,1,m),k=2,nres-2)
 259 !         write(icsa_native_int,200) (dihang_in(3,k,1,m),k=2,nres-1)
 260 !         write(icsa_native_int,200) (dihang_in(4,k,1,m),k=2,nres-1)
 261 !      enddo
 262 !     close(icsa_native_int)
 263
 264 !      do m=mm+2,n
 265 !       do i=1,4
 266 !        do j=2,nres-1
 267 !         dihang_in(i,j,1,m)=dihang_in(i,j,1,m)*deg2rad
 268 !        enddo
 269 !       enddo
 270 !      enddo
 271
 272       call dihang_to_c(dihang_in(1,1,1,1))
 273
 274 ! Store c to cref (they are in COMMON.CHAIN).
 275       do k=1,2*nres
 276        do kk=1,3
 277         crefjlee(kk,k)=c(kk,k)
 278        enddo
 279       enddo
 280
 281       call contact(.true.,ncont_ref,icont_ref,co)
 282
 283 !      do k=1,nres
 284 !       write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
 285 !      enddo
 286       close(icsa_native_int)
 287
 288   200 format (8f10.4)
 289
 290       return
 291       end subroutine from_int
 292 !-----------------------------------------------------------------------------
 293       subroutine dihang_to_c(aarray)
 294
 295       use geometry_data
 296       use csa_data
 297       use geometry, only:chainbuild
 298 !      implicit real*8 (a-h,o-z)
 299 !      include 'DIMENSIONS'
 300 !      include 'COMMON.CSA'
 301 !      include 'COMMON.BANK'
 302 !      include 'COMMON.CHAIN'
 303 !      include 'COMMON.GEO'
 304 !      include 'COMMON.VAR'
 305       integer :: i
 306       real(kind=8),dimension(mxang,nres,mxch) :: aarray !(mxang,maxres,mxch)
 307
 308 !     do i=4,nres
 309 !      phi(i)=dihang_in(1,i-2,1,1)
 310 !     enddo
 311       do i=2,nres-1
 312        theta(i+1)=aarray(1,i,1)
 313        phi(i+2)=aarray(2,i,1)
 314        alph(i)=aarray(3,i,1)
 315        omeg(i)=aarray(4,i,1)
 316       enddo
 317
 318       call chainbuild
 319
 320       return
 321       end subroutine dihang_to_c
 322 !-----------------------------------------------------------------------------
 323 ! geomout.F
 324 !-----------------------------------------------------------------------------
 325 #ifdef NOXDR
 326       subroutine cartout(time)
 327 #else
 328       subroutine cartoutx(time)
 329 #endif
 330       use geometry_data, only: c,nres
 331       use energy_data
 332       use MD_data, only: potE,t_bath
 333 !      implicit real*8 (a-h,o-z)
 334 !      include 'DIMENSIONS'
 335 !      include 'COMMON.CHAIN'
 336 !      include 'COMMON.INTERACT'
 337 !      include 'COMMON.NAMES'
 338 !      include 'COMMON.IOUNITS'
 339 !      include 'COMMON.HEADER'
 340 !      include 'COMMON.SBRIDGE'
 341 !      include 'COMMON.DISTFIT'
 342 !      include 'COMMON.MD'
 343       real(kind=8) :: time
 344 !el  local variables
 345       integer :: j,k,i
 346
 347 #if defined(AIX) || defined(PGI)
 348       open(icart,file=cartname,position="append")
 349 #else
 350       open(icart,file=cartname,access="append")
 351 #endif
 352       write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 353       if (dyn_ss) then
 354        write (icart,'(i4,$)') &
 355          nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)
 356       else
 357        write (icart,'(i4,$)') &
 358          nss,(ihpb(j),jhpb(j),j=1,nss)
 359        endif
 360        write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,&
 361        (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),&
 362        (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 363       write (icart,'(8f10.5)') &
 364        ((c(k,j),k=1,3),j=1,nres),&
 365        ((c(k,j+nres),k=1,3),j=nnt,nct)
 366       close(icart)
 367       return
 368
 369 #ifdef NOXDR
 370       end subroutine cartout
 371 #else
 372       end subroutine cartoutx
 373 #endif
 374 !-----------------------------------------------------------------------------
 375 #ifndef NOXDR
 376       subroutine cartout(time)
 377 !      implicit real*8 (a-h,o-z)
 378 !      include 'DIMENSIONS'
 379       use geometry_data, only: c,nres
 380       use energy_data
 381       use MD_data, only: potE,t_bath
 382 #ifdef MPI
 383       use MPI_data
 384       include 'mpif.h'
 385 !      include 'COMMON.SETUP'
 386 #else
 387       integer,parameter :: me=0
 388 #endif
 389 !      include 'COMMON.CHAIN'
 390 !      include 'COMMON.INTERACT'
 391 !      include 'COMMON.NAMES'
 392 !      include 'COMMON.IOUNITS'
 393 !      include 'COMMON.HEADER'
 394 !      include 'COMMON.SBRIDGE'
 395 !      include 'COMMON.DISTFIT'
 396 !      include 'COMMON.MD'
 397       real(kind=8) :: time
 398       integer :: iret,itmp
 399       real(kind=4) :: prec
 400       real(kind=4),dimension(3,2*nres+2) :: xcoord      !(3,maxres2+2)  (maxres2=2*maxres
 401 !el  local variables
 402       integer :: j,i,ixdrf
 403
 404 #ifdef AIX
 405       call xdrfopen_(ixdrf,cartname, "a", iret)
 406       call xdrffloat_(ixdrf, real(time), iret)
 407       call xdrffloat_(ixdrf, real(potE), iret)
 408       call xdrffloat_(ixdrf, real(uconst), iret)
 409       call xdrffloat_(ixdrf, real(uconst_back), iret)
 410       call xdrffloat_(ixdrf, real(t_bath), iret)
 411       call xdrfint_(ixdrf, nss, iret)
 412       do j=1,nss
 413        if (dyn_ss) then
 414         call xdrfint_(ixdrf, idssb(j)+nres, iret)
 415         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 416        else
 417         call xdrfint_(ixdrf, ihpb(j), iret)
 418         call xdrfint_(ixdrf, jhpb(j), iret)
 419        endif
 420       enddo
 421       call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
 422       do i=1,nfrag
 423         call xdrffloat_(ixdrf, real(qfrag(i)), iret)
 424       enddo
 425       do i=1,npair
 426         call xdrffloat_(ixdrf, real(qpair(i)), iret)
 427       enddo
 428       do i=1,nfrag_back
 429         call xdrffloat_(ixdrf, real(utheta(i)), iret)
 430         call xdrffloat_(ixdrf, real(ugamma(i)), iret)
 431         call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
 432       enddo
 433 #else
 434       call xdrfopen(ixdrf,cartname, "a", iret)
 435       call xdrffloat(ixdrf, real(time), iret)
 436       call xdrffloat(ixdrf, real(potE), iret)
 437       call xdrffloat(ixdrf, real(uconst), iret)
 438       call xdrffloat(ixdrf, real(uconst_back), iret)
 439       call xdrffloat(ixdrf, real(t_bath), iret)
 440       call xdrfint(ixdrf, nss, iret)
 441       do j=1,nss
 442        if (dyn_ss) then
 443         call xdrfint(ixdrf, idssb(j)+nres, iret)
 444         call xdrfint(ixdrf, jdssb(j)+nres, iret)
 445        else
 446         call xdrfint(ixdrf, ihpb(j), iret)
 447         call xdrfint(ixdrf, jhpb(j), iret)
 448        endif
 449       enddo
 450       call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
 451       do i=1,nfrag
 452         call xdrffloat(ixdrf, real(qfrag(i)), iret)
 453       enddo
 454       do i=1,npair
 455         call xdrffloat(ixdrf, real(qpair(i)), iret)
 456       enddo
 457       do i=1,nfrag_back
 458         call xdrffloat(ixdrf, real(utheta(i)), iret)
 459         call xdrffloat(ixdrf, real(ugamma(i)), iret)
 460         call xdrffloat(ixdrf, real(uscdiff(i)), iret)
 461       enddo
 462 #endif
 463       prec=10000.0
 464       do i=1,nres
 465        do j=1,3
 466         xcoord(j,i)=c(j,i)
 467        enddo
 468       enddo
 469       do i=nnt,nct
 470        do j=1,3
 471         xcoord(j,nres+i-nnt+1)=c(j,i+nres)
 472        enddo
 473       enddo
 474
 475       itmp=nres+nct-nnt+1
 476 #ifdef AIX
 477       call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
 478       call xdrfclose_(ixdrf, iret)
 479 #else
 480       call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
 481       call xdrfclose(ixdrf, iret)
 482 #endif
 483       return
 484       end subroutine cartout
 485 #endif
 486 !-----------------------------------------------------------------------------
 487       subroutine statout(itime)
 488
 489       use energy_data
 490       use control_data
 491       use MD_data
 492       use MPI_data
 493       use compare, only:rms_nac_nnc
 494 !      implicit real*8 (a-h,o-z)
 495 !      include 'DIMENSIONS'
 496 !      include 'COMMON.CONTROL'
 497 !      include 'COMMON.CHAIN'
 498 !      include 'COMMON.INTERACT'
 499 !      include 'COMMON.NAMES'
 500 !      include 'COMMON.IOUNITS'
 501 !      include 'COMMON.HEADER'
 502 !      include 'COMMON.SBRIDGE'
 503 !      include 'COMMON.DISTFIT'
 504 !      include 'COMMON.MD'
 505 !      include 'COMMON.REMD'
 506 !      include 'COMMON.SETUP'
 507       integer :: itime
 508       real(kind=8),dimension(0:n_ene) :: energia
 509 !      double precision gyrate
 510 !el      external gyrate
 511 !el      common /gucio/ cm
 512       character(len=256) :: line1,line2
 513       character(len=4) :: format1,format2
 514       character(len=30) :: format
 515 !el  local variables
 516       integer :: i,ii1,ii2,j
 517       real(kind=8) :: rms,frac,frac_nn,co,distance
 518
 519 #ifdef AIX
 520       if(itime.eq.0) then
 521        open(istat,file=statname,position="append")
 522       endif
 523 #else
 524 #ifdef PGI
 525       open(istat,file=statname,position="append")
 526 #else
 527       open(istat,file=statname,access="append")
 528 #endif
 529 #endif
 530        if (AFMlog.gt.0) then
 531        if (refstr) then
 532          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 533           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')&
 534                itime,totT,EK,potE,totE,&
 535                rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
 536                potEcomp(23),me
 537           format1="a133"
 538          else
 539 !C          print *,'A CHUJ',potEcomp(23)
 540           write (line1,'(i10,f15.2,7f12.3,i5,$)') &
 541                 itime,totT,EK,potE,totE,&
 542                 kinetic_T,t_bath,gyrate(),&
 543                 potEcomp(23),me
 544           format1="a114"
 545         endif
 546        else if (selfguide.gt.0) then
 547        distance=0.0
 548        do j=1,3
 549        distance=distance+(c(j,afmend)-c(j,afmbeg))**2
 550        enddo
 551        distance=dsqrt(distance)
 552        if (refstr) then
 553          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 554           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
 555          f9.3,i5,$)') &
 556                itime,totT,EK,potE,totE,&
 557                rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
 558                distance,potEcomp(23),me
 559           format1="a133"
 560 !C          print *,"CHUJOWO"
 561          else
 562 !C          print *,'A CHUJ',potEcomp(23)
 563           write (line1,'(i10,f15.2,8f12.3,i5,$)')&
 564                 itime,totT,EK,potE,totE, &
 565                 kinetic_T,t_bath,gyrate(),&
 566                 distance,potEcomp(23),me
 567           format1="a114"
 568         endif
 569        else
 570        if (refstr) then
 571          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 572           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
 573                 itime,totT,EK,potE,totE,&
 574                 rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
 575           format1="a133"
 576         else
 577           write (line1,'(i10,f15.2,7f12.3,i5,$)') &
 578                  itime,totT,EK,potE,totE,&
 579                  amax,kinetic_T,t_bath,gyrate(),me
 580           format1="a114"
 581         endif
 582         ENDIF
 583         if(usampl.and.totT.gt.eq_time) then
 584            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,&
 585             (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),&
 586             (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 587            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair &
 588                    +21*nfrag_back
 589         else
 590            format2="a001"
 591            line2=' '
 592         endif
 593         if (print_compon) then
 594           if(itime.eq.0) then
 595            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 596                                                            ",20a12)"
 597            write (istat,format) "#","",&
 598             (ename(print_order(i)),i=1,nprint_ene)
 599           endif
 600           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 601                                                            ",20f12.3)"
 602           write (istat,format) line1,line2,&
 603             (potEcomp(print_order(i)),i=1,nprint_ene)
 604         else
 605           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
 606           write (istat,format) line1,line2
 607         endif
 608 #if defined(AIX)
 609         call flush(istat)
 610 #else
 611         close(istat)
 612 #endif
 613       return
 614       end subroutine  statout
 615 !-----------------------------------------------------------------------------
 616 ! readrtns_CSA.F
 617 !-----------------------------------------------------------------------------
 618       subroutine readrtns
 619
 620       use control_data
 621       use energy_data
 622       use MPI_data
 623       use muca_md, only:read_muca
 624 !      implicit real*8 (a-h,o-z)
 625 !      include 'DIMENSIONS'
 626 #ifdef MPI
 627       include 'mpif.h'
 628 #endif
 629 !      include 'COMMON.SETUP'
 630 !      include 'COMMON.CONTROL'
 631 !      include 'COMMON.SBRIDGE'
 632 !      include 'COMMON.IOUNITS'
 633       logical :: file_exist
 634       integer :: i
 635 ! Read force-field parameters except weights
 636 !      call parmread
 637 ! Read job setup parameters
 638       call read_control
 639 ! Read force-field parameters except weights
 640       call parmread
 641
 642 ! Read control parameters for energy minimzation if required
 643       if (minim) call read_minim
 644 ! Read MCM control parameters if required
 645       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
 646 ! Read MD control parameters if reqjuired
 647       if (modecalc.eq.12) call read_MDpar
 648 ! Read MREMD control parameters if required
 649       if (modecalc.eq.14) then
 650          call read_MDpar
 651          call read_REMDpar
 652       endif
 653 ! Read MUCA control parameters if required
 654       if (lmuca) call read_muca
 655 ! Read CSA control parameters if required (from fort.40 if exists
 656 ! otherwise from general input file)
 657       if (modecalc.eq.8) then
 658        inquire (file="fort.40",exist=file_exist)
 659        if (.not.file_exist) call csaread
 660       endif
 661 !fmc      if (modecalc.eq.10) call mcmfread
 662 ! Read molecule information, molecule geometry, energy-term weights, and
 663 ! restraints if requested
 664       call molread
 665 ! Print restraint information
 666 #ifdef MPI
 667       if (.not. out1file .or. me.eq.king) then
 668 #endif
 669       if (nhpb.gt.nss) &
 670       write (iout,'(a,i5,a)') "The following",nhpb-nss,&
 671        " distance constraints have been imposed"
 672       do i=nss+1,nhpb
 673         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
 674       enddo
 675 #ifdef MPI
 676       endif
 677 #endif
 678 !      print *,"Processor",myrank," leaves READRTNS"
 679 !      write(iout,*) "end readrtns"
 680       return
 681       end subroutine readrtns
 682 !-----------------------------------------------------------------------------
 683       subroutine molread
 684 !
 685 ! Read molecular data.
 686 !
 687 !      use control, only: ilen
 688       use control_data
 689       use geometry_data
 690       use energy_data
 691       use energy
 692       use compare_data
 693       use MD_data, only: t_bath
 694       use MPI_data
 695       use compare, only:seq_comp,contact
 696       use control
 697 !      implicit real*8 (a-h,o-z)
 698 !      include 'DIMENSIONS'
 699 #ifdef MPI
 700       include 'mpif.h'
 701       integer :: error_msg,ierror,ierr,ierrcode
 702 #endif
 703 !      include 'COMMON.IOUNITS'
 704 !      include 'COMMON.GEO'
 705 !      include 'COMMON.VAR'
 706 !      include 'COMMON.INTERACT'
 707 !      include 'COMMON.LOCAL'
 708 !      include 'COMMON.NAMES'
 709 !      include 'COMMON.CHAIN'
 710 !      include 'COMMON.FFIELD'
 711 !      include 'COMMON.SBRIDGE'
 712 !      include 'COMMON.HEADER'
 713 !      include 'COMMON.CONTROL'
 714 !      include 'COMMON.DBASE'
 715 !      include 'COMMON.THREAD'
 716 !      include 'COMMON.CONTACTS'
 717 !      include 'COMMON.TORCNSTR'
 718 !      include 'COMMON.TIME1'
 719 !      include 'COMMON.BOUNDS'
 720 !      include 'COMMON.MD'
 721 !      include 'COMMON.SETUP'
 722       character(len=4),dimension(:,:),allocatable :: sequence   !(maxres,maxmolecules)
 723 !      integer :: rescode
 724 !      double precision x(maxvar)
 725       character(len=256) :: pdbfile
 726       character(len=800) :: weightcard
 727       character(len=80) :: weightcard_t!,ucase
 728 !      integer,dimension(:),allocatable :: itype_pdb    !(maxres)
 729 !      common /pizda/ itype_pdb
 730       logical :: fail   !seq_comp,
 731       real(kind=8) :: energia(0:n_ene)
 732 !      integer ilen
 733 !el      external ilen
 734 !el local varables
 735       integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
 736
 737       real(kind=8),dimension(3,maxres2+2) :: c_alloc
 738       real(kind=8),dimension(3,0:maxres2) :: dc_alloc
 739       integer,dimension(maxres) :: itype_alloc
 740
 741       integer :: iti,nsi,maxsi,itrial,itmp
 742       real(kind=8) :: wlong,scalscp,co,ssscale
 743       allocate(weights(n_ene))
 744 !-----------------------------
 745       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
 746       allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
 747       allocate(itype(maxres,5)) !(maxres)
 748       allocate(istype(maxres))
 749 !
 750 ! Zero out tables.
 751 !
 752       c(:,:)=0.0D0
 753       dc(:,:)=0.0D0
 754       itype(:,:)=0
 755 !-----------------------------
 756 !
 757 ! Body
 758 !
 759 ! Read weights of the subsequent energy terms.
 760       call card_concat(weightcard,.true.)
 761       call reada(weightcard,'WLONG',wlong,1.0D0)
 762       call reada(weightcard,'WSC',wsc,wlong)
 763       call reada(weightcard,'WSCP',wscp,wlong)
 764       call reada(weightcard,'WELEC',welec,1.0D0)
 765       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 766       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 767       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 768       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 769       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 770       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 771       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 772       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 773       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 774       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 775       call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
 776       call reada(weightcard,'WELPP',welpp,0.0d0)
 777       call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
 778       call reada(weightcard,'WELPSB',welpsb,0.0D0)
 779       call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
 780       call reada(weightcard,'WELSB',welsb,0.0D0)
 781       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
 782       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
 783       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
 784       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
 785       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
 786       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
 787       call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,0.0D0)
 788       call reada(weightcard,'WBOND',wbond,1.0D0)
 789       call reada(weightcard,'WTOR',wtor,1.0D0)
 790       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 791       call reada(weightcard,'WSHIELD',wshield,0.05D0)
 792       call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
 793       call reada(weightcard,'WLT',wliptran,1.0D0)
 794       call reada(weightcard,'WTUBE',wtube,1.0d0)
 795       call reada(weightcard,'WANG',wang,1.0D0)
 796       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 797       call reada(weightcard,'SCAL14',scal14,0.4D0)
 798       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 799       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 800       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 801       call reada(weightcard,'TEMP0',temp0,300.0d0)
 802       call reada(weightcard,'WSCBASE',wscbase,0.0D0)
 803       if (index(weightcard,'SOFT').gt.0) ipot=6
 804       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
 805       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
 806       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
 807       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
 808       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
 809       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
 810
 811 ! 12/1/95 Added weight for the multi-body term WCORR
 812       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 813       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 814       weights(1)=wsc
 815       weights(2)=wscp
 816       weights(3)=welec
 817       weights(4)=wcorr
 818       weights(5)=wcorr5
 819       weights(6)=wcorr6
 820       weights(7)=wel_loc
 821       weights(8)=wturn3
 822       weights(9)=wturn4
 823       weights(10)=wturn6
 824       weights(11)=wang
 825       weights(12)=wscloc
 826       weights(13)=wtor
 827       weights(14)=wtor_d
 828       weights(15)=wstrain
 829       weights(16)=wvdwpp
 830       weights(17)=wbond
 831       weights(18)=scal14
 832       weights(21)=wsccor
 833       if(me.eq.king.or..not.out1file) &
 834        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 835         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 836         wturn4,wturn6
 837    10 format (/'Energy-term weights (unscaled):'// &
 838        'WSCC=   ',f10.6,' (SC-SC)'/ &
 839        'WSCP=   ',f10.6,' (SC-p)'/ &
 840        'WELEC=  ',f10.6,' (p-p electr)'/ &
 841        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 842        'WBOND=  ',f10.6,' (stretching)'/ &
 843        'WANG=   ',f10.6,' (bending)'/ &
 844        'WSCLOC= ',f10.6,' (SC local)'/ &
 845        'WTOR=   ',f10.6,' (torsional)'/ &
 846        'WTORD=  ',f10.6,' (double torsional)'/ &
 847        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 848        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 849        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 850        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 851        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 852        'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
 853        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 854        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 855        'WTURN6= ',f10.6,' (turns, 6th order)')
 856       if(me.eq.king.or..not.out1file)then
 857        if (wcorr4.gt.0.0d0) then
 858         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
 859          'between contact pairs of peptide groups'
 860         write (iout,'(2(a,f5.3/))') &
 861         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
 862         'Range of quenching the correlation terms:',2*delt_corr
 863        else if (wcorr.gt.0.0d0) then
 864         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
 865          'between contact pairs of peptide groups'
 866        endif
 867        write (iout,'(a,f8.3)') &
 868         'Scaling factor of 1,4 SC-p interactions:',scal14
 869        write (iout,'(a,f8.3)') &
 870         'General scaling factor of SC-p interactions:',scalscp
 871       endif
 872       r0_corr=cutoff_corr-delt_corr
 873       do i=1,ntyp
 874         aad(i,1)=scalscp*aad(i,1)
 875         aad(i,2)=scalscp*aad(i,2)
 876         bad(i,1)=scalscp*bad(i,1)
 877         bad(i,2)=scalscp*bad(i,2)
 878       enddo
 879       call rescale_weights(t_bath)
 880       if(me.eq.king.or..not.out1file) &
 881        write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 882         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 883         wturn4,wturn6
 884    22 format (/'Energy-term weights (scaled):'// &
 885        'WSCC=   ',f10.6,' (SC-SC)'/ &
 886        'WSCP=   ',f10.6,' (SC-p)'/ &
 887        'WELEC=  ',f10.6,' (p-p electr)'/ &
 888        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 889        'WBOND=  ',f10.6,' (stretching)'/ &
 890        'WANG=   ',f10.6,' (bending)'/ &
 891        'WSCLOC= ',f10.6,' (SC local)'/ &
 892        'WTOR=   ',f10.6,' (torsional)'/ &
 893        'WTORD=  ',f10.6,' (double torsional)'/ &
 894        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 895        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 896        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 897        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 898        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 899        'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
 900        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 901        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 902        'WTURN6= ',f10.6,' (turns, 6th order)')
 903       if(me.eq.king.or..not.out1file) &
 904        write (iout,*) "Reference temperature for weights calculation:",&
 905         temp0
 906       call reada(weightcard,"D0CM",d0cm,3.78d0)
 907       call reada(weightcard,"AKCM",akcm,15.1d0)
 908       call reada(weightcard,"AKTH",akth,11.0d0)
 909       call reada(weightcard,"AKCT",akct,12.0d0)
 910       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 911       call reada(weightcard,"V2SS",v2ss,7.61d0)
 912       call reada(weightcard,"V3SS",v3ss,13.7d0)
 913       call reada(weightcard,"EBR",ebr,-5.50D0)
 914       call reada(weightcard,"ATRISS",atriss,0.301D0)
 915       call reada(weightcard,"BTRISS",btriss,0.021D0)
 916       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 917       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 918       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
 919       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 920
 921       call reada(weightcard,"HT",Ht,0.0D0)
 922       if (dyn_ss) then
 923        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
 924         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
 925         akcm=akcm/wsc*ssscale
 926         akth=akth/wsc*ssscale
 927         akct=akct/wsc*ssscale
 928         v1ss=v1ss/wsc*ssscale
 929         v2ss=v2ss/wsc*ssscale
 930         v3ss=v3ss/wsc*ssscale
 931       else
 932         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 933       endif
 934
 935       if(me.eq.king.or..not.out1file) then
 936        write (iout,*) "Parameters of the SS-bond potential:"
 937        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
 938        " AKCT",akct
 939        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 940        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 941        write (iout,*)" HT",Ht
 942        print *,'indpdb=',indpdb,' pdbref=',pdbref
 943       endif
 944       if (indpdb.gt.0 .or. pdbref) then
 945         read(inp,'(a)') pdbfile
 946         if(me.eq.king.or..not.out1file) &
 947          write (iout,'(2a)') 'PDB data will be read from file ',&
 948          pdbfile(:ilen(pdbfile))
 949         open(ipdbin,file=pdbfile,status='old',err=33)
 950         goto 34
 951   33    write (iout,'(a)') 'Error opening PDB file.'
 952         stop
 953   34    continue
 954 !        print *,'Begin reading pdb data'
 955         call readpdb
 956 !        print *,'Finished reading pdb data'
 957         if(me.eq.king.or..not.out1file) &
 958          write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
 959          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
 960 !el        if(.not.allocated(itype_pdb))
 961         allocate(itype_pdb(nres))
 962         do i=1,nres
 963           itype_pdb(i)=itype(i,1)
 964         enddo
 965         close (ipdbin)
 966         nnt=nstart_sup
 967         nct=nstart_sup+nsup-1
 968 !el        if(.not.allocated(icont_ref))
 969         allocate(icont_ref(2,(nres/2)*nres)) ! maxcont=12*maxres
 970         call contact(.false.,ncont_ref,icont_ref,co)
 971
 972         if (sideadd) then
 973          if(me.eq.king.or..not.out1file) &
 974           write(iout,*)'Adding sidechains'
 975          maxsi=1000
 976          do i=2,nres-1
 977           iti=itype(i,1)
 978           if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
 979             nsi=0
 980             fail=.true.
 981             do while (fail.and.nsi.le.maxsi)
 982               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 983               nsi=nsi+1
 984             enddo
 985             if(fail) write(iout,*)'Adding sidechain failed for res ',&
 986                     i,' after ',nsi,' trials'
 987           endif
 988          enddo
 989         endif
 990       endif
 991
 992       if (indpdb.eq.0) then
 993       nres_molec(:)=0
 994         allocate(sequence(maxres,5))
 995 !      itype(:,:)=0
 996       itmp=0
 997       if (protein) then
 998 ! Read sequence if not taken from the pdb file.
 999         molec=1
1000         read (inp,*) nres_molec(molec)
1001         print *,'nres=',nres
1002         if (iscode.gt.0) then
1003           read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
1004         else
1005           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1006         endif
1007 !        read(inp,*) weightcard_t
1008 !        print *,"po seq" weightcard_t
1009 ! Convert sequence to numeric code
1010
1011         do i=1,nres_molec(molec)
1012           itmp=itmp+1
1013           itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1014           print *,itype(i,1)
1015
1016         enddo
1017        endif
1018 !        read(inp,*) weightcard_t
1019 !        print *,"po seq", weightcard_t
1020
1021        if (nucleic) then
1022 ! Read sequence if not taken from the pdb file.
1023         molec=2
1024         read (inp,*) nres_molec(molec)
1025 !        print *,'nres=',nres
1026 !        allocate(sequence(maxres,5))
1027 !        if (iscode.gt.0) then
1028           read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
1029 ! Convert sequence to numeric code
1030
1031         do i=1,nres_molec(molec)
1032           itmp=itmp+1
1033           istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)
1034           itype(itmp,molec)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
1035         enddo
1036        endif
1037
1038        if (ions) then
1039 ! Read sequence if not taken from the pdb file.
1040         molec=5
1041         read (inp,*) nres_molec(molec)
1042 !        print *,'nres=',nres
1043           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1044 ! Convert sequence to numeric code
1045         print *,nres_molec(molec)
1046         do i=1,nres_molec(molec)
1047           itmp=itmp+1
1048           print *,itmp,"itmp"
1049           itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1050         enddo
1051        endif
1052        nres=0
1053        do i=1,5
1054         nres=nres+nres_molec(i)
1055         print *,nres_molec(i)
1056        enddo
1057
1058 ! Assign initial virtual bond lengths
1059         if(.not.allocated(molnum)) then
1060          allocate(molnum(nres+1))
1061          itmp=0
1062         do i=1,5
1063                do j=1,nres_molec(i)
1064                itmp=itmp+1
1065               molnum(itmp)=i
1066                enddo
1067          enddo
1068 !        print *,nres_molec(i)
1069         endif
1070         if(.not.allocated(vbld)) allocate(vbld(2*nres))
1071         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1072         do i=2,nres
1073           vbld(i)=vbl
1074           vbld_inv(i)=vblinv
1075         enddo
1076         do i=2,nres-1
1077           print *, "molnum",molnum(i),itype(i,molnum(i))
1078           vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
1079           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
1080 !          write (iout,*) "i",i," itype",itype(i,1),
1081 !     &      " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1082         enddo
1083       endif
1084 !      print *,nres
1085 !      print '(20i4)',(itype(i,1),i=1,nres)
1086 !----------------------------
1087 !el reallocate tables
1088 !      do i=1,maxres2
1089 !        do j=1,3
1090 !          c_alloc(j,i)=c(j,i)
1091 !          dc_alloc(j,i)=dc(j,i)
1092 !        enddo
1093 !      enddo
1094 !      do i=1,maxres
1095 !elwrite(iout,*) "itype",i,itype(i,1)
1096 !        itype_alloc(i)=itype(i,1)
1097 !      enddo
1098
1099 !      deallocate(c)
1100 !      deallocate(dc)
1101 !      deallocate(itype)
1102 !      allocate(c(3,2*nres+4))
1103 !      allocate(dc(3,0:2*nres+2))
1104 !      allocate(itype(nres+2))
1105       allocate(itel(nres+2))
1106       itel(:)=0
1107
1108 !      do i=1,2*nres+2
1109 !        do j=1,3
1110 !          c(j,i)=c_alloc(j,i)
1111 !          dc(j,i)=dc_alloc(j,i)
1112 !        enddo
1113 !      enddo
1114 !      do i=1,nres+2
1115 !        itype(i,1)=itype_alloc(i)
1116 !        itel(i)=0
1117 !      enddo
1118 !--------------------------
1119       do i=1,nres
1120 #ifdef PROCOR
1121         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1122 #else
1123         if (itype(i,1).eq.ntyp1) then
1124 #endif
1125           itel(i)=0
1126 #ifdef PROCOR
1127         else if (iabs(itype(i+1,1)).ne.20) then
1128 #else
1129         else if (iabs(itype(i,1)).ne.20) then
1130 #endif
1131           itel(i)=1
1132         else
1133           itel(i)=2
1134         endif
1135       enddo
1136       if(me.eq.king.or..not.out1file)then
1137        write (iout,*) "ITEL"
1138        print *,nres,"nres"
1139        do i=1,nres-1
1140          write (iout,*) i,itype(i,1),itel(i)
1141        enddo
1142        print *,'Call Read_Bridge.'
1143       endif
1144       call read_bridge
1145 !--------------------------------
1146        print *,"tu dochodze"
1147 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1148       call alloc_geo_arrays
1149       call alloc_ener_arrays
1150 !--------------------------------
1151 ! 8/13/98 Set limits to generating the dihedral angles
1152       do i=1,nres
1153         phibound(1,i)=-pi
1154         phibound(2,i)=pi
1155       enddo
1156       read (inp,*) ndih_constr
1157       if (ndih_constr.gt.0) then
1158         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1159         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1160         allocate(ftors(ndih_constr)) !(maxdih_constr)
1161
1162 !        read (inp,*) ftors
1163         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), &
1164         i=1,ndih_constr)
1165         if(me.eq.king.or..not.out1file)then
1166          write (iout,*) &
1167          'There are',ndih_constr,' constraints on phi angles.'
1168          do i=1,ndih_constr
1169           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), &
1170           ftors(i)
1171          enddo
1172         endif
1173         do i=1,ndih_constr
1174           phi0(i)=deg2rad*phi0(i)
1175           drange(i)=deg2rad*drange(i)
1176         enddo
1177 !        if(me.eq.king.or..not.out1file) &
1178 !         write (iout,*) 'FTORS',ftors
1179         do i=1,ndih_constr
1180           ii = idih_constr(i)
1181           phibound(1,ii) = phi0(i)-drange(i)
1182           phibound(2,ii) = phi0(i)+drange(i)
1183         enddo
1184       endif
1185       if (with_theta_constr) then
1186 !C with_theta_constr is keyword allowing for occurance of theta constrains
1187       read (inp,*) ntheta_constr
1188 !C ntheta_constr is the number of theta constrains
1189       if (ntheta_constr.gt.0) then
1190 !C        read (inp,*) ftors
1191         allocate(itheta_constr(ntheta_constr))
1192         allocate(theta_constr0(ntheta_constr))
1193         allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1194         read (inp,*) (itheta_constr(i),theta_constr0(i), &
1195        theta_drange(i),for_thet_constr(i), &
1196        i=1,ntheta_constr)
1197 !C the above code reads from 1 to ntheta_constr
1198 !C itheta_constr(i) residue i for which is theta_constr
1199 !C theta_constr0 the global minimum value
1200 !C theta_drange is range for which there is no energy penalty
1201 !C for_thet_constr is the force constant for quartic energy penalty
1202 !C E=k*x**4
1203         if(me.eq.king.or..not.out1file)then
1204          write (iout,*) &
1205         'There are',ntheta_constr,' constraints on phi angles.'
1206          do i=1,ntheta_constr
1207           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1208          theta_drange(i), &
1209          for_thet_constr(i)
1210          enddo
1211         endif
1212         do i=1,ntheta_constr
1213           theta_constr0(i)=deg2rad*theta_constr0(i)
1214           theta_drange(i)=deg2rad*theta_drange(i)
1215         enddo
1216 !C        if(me.eq.king.or..not.out1file)
1217 !C     &   write (iout,*) 'FTORS',ftors
1218 !C        do i=1,ntheta_constr
1219 !C          ii = itheta_constr(i)
1220 !C          thetabound(1,ii) = phi0(i)-drange(i)
1221 !C          thetabound(2,ii) = phi0(i)+drange(i)
1222 !C        enddo
1223       endif ! ntheta_constr.gt.0
1224       endif! with_theta_constr
1225
1226       nnt=1
1227 #ifdef MPI
1228       if (me.eq.king) then
1229 #endif
1230        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1231        do i=1,nres
1232          write (iout,'(a3,i5,2f10.1)') &
1233          restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1234        enddo
1235 #ifdef MP
1236       endif
1237 #endif
1238       nct=nres
1239       print *,'NNT=',NNT,' NCT=',NCT
1240       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1241       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1242       if (pdbref) then
1243         if(me.eq.king.or..not.out1file) &
1244          write (iout,'(a,i3)') 'nsup=',nsup
1245         nstart_seq=nnt
1246         if (nsup.le.(nct-nnt+1)) then
1247           do i=0,nct-nnt+1-nsup
1248             if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1249               nstart_seq=nnt+i
1250               goto 111
1251             endif
1252           enddo
1253           write (iout,'(a)') &
1254                   'Error - sequences to be superposed do not match.'
1255           stop
1256         else
1257           do i=0,nsup-(nct-nnt+1)
1258             if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1259             then
1260               nstart_sup=nstart_sup+i
1261               nsup=nct-nnt+1
1262               goto 111
1263             endif
1264           enddo
1265           write (iout,'(a)') &
1266                   'Error - sequences to be superposed do not match.'
1267         endif
1268   111   continue
1269         if (nsup.eq.0) nsup=nct-nnt
1270         if (nstart_sup.eq.0) nstart_sup=nnt
1271         if (nstart_seq.eq.0) nstart_seq=nnt
1272         if(me.eq.king.or..not.out1file) &
1273          write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1274                        ' nstart_seq=',nstart_seq !,"242343453254"
1275       endif
1276 !--- Zscore rms -------
1277       if (nz_start.eq.0) nz_start=nnt
1278       if (nz_end.eq.0 .and. nsup.gt.0) then
1279         nz_end=nnt+nsup-1
1280       else if (nz_end.eq.0) then
1281         nz_end=nct
1282       endif
1283       if(me.eq.king.or..not.out1file)then
1284        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1285        write (iout,*) 'IZ_SC=',iz_sc
1286       endif
1287 !----------------------
1288       call init_int_table
1289       if (refstr) then
1290         if (.not.pdbref) then
1291           call read_angles(inp,*38)
1292           goto 39
1293    38     write (iout,'(a)') 'Error reading reference structure.'
1294 #ifdef MPI
1295           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1296           stop 'Error reading reference structure'
1297 #endif
1298    39     call chainbuild
1299           call setup_var
1300 !zscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1301           nstart_sup=nnt
1302           nstart_seq=nnt
1303           nsup=nct-nnt+1
1304           kkk=1
1305           do i=1,2*nres
1306             do j=1,3
1307               cref(j,i,kkk)=c(j,i)
1308             enddo
1309           enddo
1310           call contact(.true.,ncont_ref,icont_ref,co)
1311         endif
1312 !        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1313 !        call flush(iout)
1314 !EL        if (constr_dist.gt.0) call read_dist_constr
1315 !EL        write (iout,*) "After read_dist_constr nhpb",nhpb
1316 !EL        if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1317 !EL        call hpb_partition
1318         if(me.eq.king.or..not.out1file) &
1319          write (iout,*) 'Contact order:',co
1320         if (pdbref) then
1321         if(me.eq.king.or..not.out1file) &
1322          write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1323         do i=1,ncont_ref
1324           do j=1,2
1325             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1326           enddo
1327           if(me.eq.king.or..not.out1file) &
1328            write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1329            icont_ref(1,i),' ',&
1330            restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1331         enddo
1332         endif
1333       endif
1334         if (constr_dist.gt.0) call read_dist_constr
1335         write (iout,*) "After read_dist_constr nhpb",nhpb
1336         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1337         call hpb_partition
1338
1339       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1340           .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1341           modecalc.ne.10) then
1342 ! If input structure hasn't been supplied from the PDB file read or generate
1343 ! initial geometry.
1344         if (iranconf.eq.0 .and. .not. extconf) then
1345           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1346            write (iout,'(a)') 'Initial geometry will be read in.'
1347           if (read_cart) then
1348             read(inp,'(8f10.5)',end=36,err=36) &
1349              ((c(l,k),l=1,3),k=1,nres),&
1350              ((c(l,k+nres),l=1,3),k=nnt,nct)
1351             write (iout,*) "Exit READ_CART"
1352             write (iout,'(8f10.5)') &
1353              ((c(l,k),l=1,3),k=1,nres)
1354             write (iout,'(8f10.5)') &
1355              ((c(l,k+nres),l=1,3),k=nnt,nct)
1356             call int_from_cart1(.true.)
1357             write (iout,*) "Finish INT_TO_CART"
1358             do i=1,nres-1
1359               do j=1,3
1360                 dc(j,i)=c(j,i+1)-c(j,i)
1361                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1362               enddo
1363             enddo
1364             do i=nnt,nct
1365               if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1366                 do j=1,3
1367                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1368                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1369                 enddo
1370               endif
1371             enddo
1372             return
1373           else
1374             call read_angles(inp,*36)
1375           endif
1376           goto 37
1377    36     write (iout,'(a)') 'Error reading angle file.'
1378 #ifdef MPI
1379           call mpi_finalize( MPI_COMM_WORLD,IERR )
1380 #endif
1381           stop 'Error reading angle file.'
1382    37     continue
1383         else if (extconf) then
1384          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1385           write (iout,'(a)') 'Extended chain initial geometry.'
1386          do i=3,nres
1387           theta(i)=90d0*deg2rad
1388          enddo
1389          do i=4,nres
1390           phi(i)=180d0*deg2rad
1391          enddo
1392          do i=2,nres-1
1393           alph(i)=110d0*deg2rad
1394          enddo
1395          do i=2,nres-1
1396           omeg(i)=-120d0*deg2rad
1397           if (itype(i,1).le.0) omeg(i)=-omeg(i)
1398          enddo
1399          call chainbuild
1400         else
1401           if(me.eq.king.or..not.out1file) &
1402            write (iout,'(a)') 'Random-generated initial geometry.'
1403
1404
1405 #ifdef MPI
1406           if (me.eq.king  .or. fg_rank.eq.0 .and. &
1407                  ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1408 #endif
1409             do itrial=1,100
1410               itmp=1
1411               call gen_rand_conf(itmp,*30)
1412               goto 40
1413    30         write (iout,*) 'Failed to generate random conformation',&
1414                 ', itrial=',itrial
1415               write (*,*) 'Processor:',me,&
1416                 ' Failed to generate random conformation',&
1417                 ' itrial=',itrial
1418               call intout
1419
1420 #ifdef AIX
1421               call flush_(iout)
1422 #else
1423               call flush(iout)
1424 #endif
1425             enddo
1426             write (iout,'(a,i3,a)') 'Processor:',me,&
1427               ' error in generating random conformation.'
1428             write (*,'(a,i3,a)') 'Processor:',me,&
1429               ' error in generating random conformation.'
1430             call flush(iout)
1431 #ifdef MPI
1432             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1433    40       continue
1434           endif
1435 #else
1436           do itrial=1,100
1437             itmp=1
1438             call gen_rand_conf(itmp,*335)
1439             goto 40
1440   335       write (iout,*) 'Failed to generate random conformation',&
1441               ', itrial=',itrial
1442             write (*,*) 'Failed to generate random conformation',&
1443               ', itrial=',itrial
1444           enddo
1445           write (iout,'(a,i3,a)') 'Processor:',me,&
1446             ' error in generating random conformation.'
1447           write (*,'(a,i3,a)') 'Processor:',me,&
1448             ' error in generating random conformation.'
1449           stop
1450    40     continue
1451 #endif
1452         endif
1453       elseif (modecalc.eq.4) then
1454         read (inp,'(a)') intinname
1455         open (intin,file=intinname,status='old',err=333)
1456         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1457         write (iout,'(a)') 'intinname',intinname
1458         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1459         goto 334
1460   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1461 #ifdef MPI
1462         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1463 #endif
1464         stop 'Error opening angle file.'
1465   334   continue
1466
1467       endif
1468 ! Generate distance constraints, if the PDB structure is to be regularized.
1469       if (nthread.gt.0) then
1470         call read_threadbase
1471       endif
1472       call setup_var
1473       if (me.eq.king .or. .not. out1file) &
1474        call intout
1475       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1476         write (iout,'(/a,i3,a)') &
1477         'The chain contains',ns,' disulfide-bridging cysteines.'
1478         write (iout,'(20i4)') (iss(i),i=1,ns)
1479        if (dyn_ss) then
1480           write(iout,*)"Running with dynamic disulfide-bond formation"
1481        else
1482         write (iout,'(/a/)') 'Pre-formed links are:'
1483         do i=1,nss
1484           i1=ihpb(i)-nres
1485           i2=jhpb(i)-nres
1486           it1=itype(i1,1)
1487           it2=itype(i2,1)
1488           if (me.eq.king.or..not.out1file) &
1489           write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1490           restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1491           ebr,forcon(i)
1492         enddo
1493         write (iout,'(a)')
1494        endif
1495       endif
1496       if (ns.gt.0.and.dyn_ss) then
1497           do i=nss+1,nhpb
1498             ihpb(i-nss)=ihpb(i)
1499             jhpb(i-nss)=jhpb(i)
1500             forcon(i-nss)=forcon(i)
1501             dhpb(i-nss)=dhpb(i)
1502           enddo
1503           nhpb=nhpb-nss
1504           nss=0
1505           call hpb_partition
1506           do i=1,ns
1507             dyn_ss_mask(iss(i))=.true.
1508           enddo
1509       endif
1510       if (i2ndstr.gt.0) call secstrp2dihc
1511 !      call geom_to_var(nvar,x)
1512 !      call etotal(energia(0))
1513 !      call enerprint(energia(0))
1514 !      call briefout(0,etot)
1515 !      stop
1516 !d    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1517 !d    write (iout,'(a)') 'Variable list:'
1518 !d    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1519 #ifdef MPI
1520       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1521         write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1522         'Processor',myrank,': end reading molecular data.'
1523 #endif
1524       return
1525       end subroutine molread
1526 !-----------------------------------------------------------------------------
1527 !-----------------------------------------------------------------------------
1528       end module io