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 control parameters for energy minimzation if required
 640       if (minim) call read_minim
 641 ! Read MCM control parameters if required
 642       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
 643 ! Read MD control parameters if reqjuired
 644       if (modecalc.eq.12) call read_MDpar
 645 ! Read MREMD control parameters if required
 646       if (modecalc.eq.14) then
 647          call read_MDpar
 648          call read_REMDpar
 649       endif
 650 ! Read MUCA control parameters if required
 651       if (lmuca) call read_muca
 652 ! Read CSA control parameters if required (from fort.40 if exists
 653 ! otherwise from general input file)
 654       if (modecalc.eq.8) then
 655        inquire (file="fort.40",exist=file_exist)
 656        if (.not.file_exist) call csaread
 657       endif
 658 !fmc      if (modecalc.eq.10) call mcmfread
 659 ! Read molecule information, molecule geometry, energy-term weights, and
 660 ! restraints if requested
 661       call molread
 662 ! Print restraint information
 663 #ifdef MPI
 664       if (.not. out1file .or. me.eq.king) then
 665 #endif
 666       if (nhpb.gt.nss) &
 667       write (iout,'(a,i5,a)') "The following",nhpb-nss,&
 668        " distance constraints have been imposed"
 669       do i=nss+1,nhpb
 670         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
 671       enddo
 672 #ifdef MPI
 673       endif
 674 #endif
 675 !      print *,"Processor",myrank," leaves READRTNS"
 676 !      write(iout,*) "end readrtns"
 677       return
 678       end subroutine readrtns
 679 !-----------------------------------------------------------------------------
 680       subroutine molread
 681 !
 682 ! Read molecular data.
 683 !
 684 !      use control, only: ilen
 685       use control_data
 686       use geometry_data
 687       use energy_data
 688       use energy
 689       use compare_data
 690       use MD_data, only: t_bath
 691       use MPI_data
 692       use compare, only:seq_comp,contact
 693       use control
 694 !      implicit real*8 (a-h,o-z)
 695 !      include 'DIMENSIONS'
 696 #ifdef MPI
 697       include 'mpif.h'
 698       integer :: error_msg,ierror,ierr,ierrcode
 699 #endif
 700 !      include 'COMMON.IOUNITS'
 701 !      include 'COMMON.GEO'
 702 !      include 'COMMON.VAR'
 703 !      include 'COMMON.INTERACT'
 704 !      include 'COMMON.LOCAL'
 705 !      include 'COMMON.NAMES'
 706 !      include 'COMMON.CHAIN'
 707 !      include 'COMMON.FFIELD'
 708 !      include 'COMMON.SBRIDGE'
 709 !      include 'COMMON.HEADER'
 710 !      include 'COMMON.CONTROL'
 711 !      include 'COMMON.DBASE'
 712 !      include 'COMMON.THREAD'
 713 !      include 'COMMON.CONTACTS'
 714 !      include 'COMMON.TORCNSTR'
 715 !      include 'COMMON.TIME1'
 716 !      include 'COMMON.BOUNDS'
 717 !      include 'COMMON.MD'
 718 !      include 'COMMON.SETUP'
 719       character(len=4),dimension(:,:),allocatable :: sequence   !(maxres,maxmolecules)
 720 !      integer :: rescode
 721 !      double precision x(maxvar)
 722       character(len=256) :: pdbfile
 723       character(len=480) :: weightcard
 724       character(len=80) :: weightcard_t!,ucase
 725 !      integer,dimension(:),allocatable :: itype_pdb    !(maxres)
 726 !      common /pizda/ itype_pdb
 727       logical :: fail   !seq_comp,
 728       real(kind=8) :: energia(0:n_ene)
 729 !      integer ilen
 730 !el      external ilen
 731 !el local varables
 732       integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
 733
 734       real(kind=8),dimension(3,maxres2+2) :: c_alloc
 735       real(kind=8),dimension(3,0:maxres2) :: dc_alloc
 736       integer,dimension(maxres) :: itype_alloc
 737
 738       integer :: iti,nsi,maxsi,itrial,itmp
 739       real(kind=8) :: wlong,scalscp,co,ssscale
 740       allocate(weights(n_ene))
 741 !-----------------------------
 742       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
 743       allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
 744       allocate(itype(maxres,5)) !(maxres)
 745       allocate(istype(maxres))
 746 !
 747 ! Zero out tables.
 748 !
 749       c(:,:)=0.0D0
 750       dc(:,:)=0.0D0
 751       itype(:,:)=0
 752 !-----------------------------
 753 !
 754 ! Body
 755 !
 756 ! Read weights of the subsequent energy terms.
 757       call card_concat(weightcard,.true.)
 758       call reada(weightcard,'WLONG',wlong,1.0D0)
 759       call reada(weightcard,'WSC',wsc,wlong)
 760       call reada(weightcard,'WSCP',wscp,wlong)
 761       call reada(weightcard,'WELEC',welec,1.0D0)
 762       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 763       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 764       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 765       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 766       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 767       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 768       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 769       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 770       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 771       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 772       call reada(weightcard,'WBOND',wbond,1.0D0)
 773       call reada(weightcard,'WTOR',wtor,1.0D0)
 774       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 775       call reada(weightcard,'WSHIELD',wshield,0.05D0)
 776       call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
 777       call reada(weightcard,'WLT',wliptran,1.0D0)
 778       call reada(weightcard,'WTUBE',wtube,1.0d0)
 779       call reada(weightcard,'WANG',wang,1.0D0)
 780       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 781       call reada(weightcard,'SCAL14',scal14,0.4D0)
 782       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 783       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 784       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 785       call reada(weightcard,'TEMP0',temp0,300.0d0)
 786       if (index(weightcard,'SOFT').gt.0) ipot=6
 787       call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0)
 788 ! 12/1/95 Added weight for the multi-body term WCORR
 789       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 790       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 791       weights(1)=wsc
 792       weights(2)=wscp
 793       weights(3)=welec
 794       weights(4)=wcorr
 795       weights(5)=wcorr5
 796       weights(6)=wcorr6
 797       weights(7)=wel_loc
 798       weights(8)=wturn3
 799       weights(9)=wturn4
 800       weights(10)=wturn6
 801       weights(11)=wang
 802       weights(12)=wscloc
 803       weights(13)=wtor
 804       weights(14)=wtor_d
 805       weights(15)=wstrain
 806       weights(16)=wvdwpp
 807       weights(17)=wbond
 808       weights(18)=scal14
 809       weights(21)=wsccor
 810       if(me.eq.king.or..not.out1file) &
 811        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 812         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 813         wturn4,wturn6
 814    10 format (/'Energy-term weights (unscaled):'// &
 815        'WSCC=   ',f10.6,' (SC-SC)'/ &
 816        'WSCP=   ',f10.6,' (SC-p)'/ &
 817        'WELEC=  ',f10.6,' (p-p electr)'/ &
 818        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 819        'WBOND=  ',f10.6,' (stretching)'/ &
 820        'WANG=   ',f10.6,' (bending)'/ &
 821        'WSCLOC= ',f10.6,' (SC local)'/ &
 822        'WTOR=   ',f10.6,' (torsional)'/ &
 823        'WTORD=  ',f10.6,' (double torsional)'/ &
 824        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 825        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 826        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 827        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 828        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 829        'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
 830        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 831        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 832        'WTURN6= ',f10.6,' (turns, 6th order)')
 833       if(me.eq.king.or..not.out1file)then
 834        if (wcorr4.gt.0.0d0) then
 835         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
 836          'between contact pairs of peptide groups'
 837         write (iout,'(2(a,f5.3/))') &
 838         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
 839         'Range of quenching the correlation terms:',2*delt_corr
 840        else if (wcorr.gt.0.0d0) then
 841         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
 842          'between contact pairs of peptide groups'
 843        endif
 844        write (iout,'(a,f8.3)') &
 845         'Scaling factor of 1,4 SC-p interactions:',scal14
 846        write (iout,'(a,f8.3)') &
 847         'General scaling factor of SC-p interactions:',scalscp
 848       endif
 849       r0_corr=cutoff_corr-delt_corr
 850       do i=1,ntyp
 851         aad(i,1)=scalscp*aad(i,1)
 852         aad(i,2)=scalscp*aad(i,2)
 853         bad(i,1)=scalscp*bad(i,1)
 854         bad(i,2)=scalscp*bad(i,2)
 855       enddo
 856       call rescale_weights(t_bath)
 857       if(me.eq.king.or..not.out1file) &
 858        write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 859         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 860         wturn4,wturn6
 861    22 format (/'Energy-term weights (scaled):'// &
 862        'WSCC=   ',f10.6,' (SC-SC)'/ &
 863        'WSCP=   ',f10.6,' (SC-p)'/ &
 864        'WELEC=  ',f10.6,' (p-p electr)'/ &
 865        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 866        'WBOND=  ',f10.6,' (stretching)'/ &
 867        'WANG=   ',f10.6,' (bending)'/ &
 868        'WSCLOC= ',f10.6,' (SC local)'/ &
 869        'WTOR=   ',f10.6,' (torsional)'/ &
 870        'WTORD=  ',f10.6,' (double torsional)'/ &
 871        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 872        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 873        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 874        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 875        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 876        'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
 877        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 878        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 879        'WTURN6= ',f10.6,' (turns, 6th order)')
 880       if(me.eq.king.or..not.out1file) &
 881        write (iout,*) "Reference temperature for weights calculation:",&
 882         temp0
 883       call reada(weightcard,"D0CM",d0cm,3.78d0)
 884       call reada(weightcard,"AKCM",akcm,15.1d0)
 885       call reada(weightcard,"AKTH",akth,11.0d0)
 886       call reada(weightcard,"AKCT",akct,12.0d0)
 887       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 888       call reada(weightcard,"V2SS",v2ss,7.61d0)
 889       call reada(weightcard,"V3SS",v3ss,13.7d0)
 890       call reada(weightcard,"EBR",ebr,-5.50D0)
 891       call reada(weightcard,"ATRISS",atriss,0.301D0)
 892       call reada(weightcard,"BTRISS",btriss,0.021D0)
 893       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 894       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 895       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
 896       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 897
 898       call reada(weightcard,"HT",Ht,0.0D0)
 899       if (dyn_ss) then
 900        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
 901         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
 902         akcm=akcm/wsc*ssscale
 903         akth=akth/wsc*ssscale
 904         akct=akct/wsc*ssscale
 905         v1ss=v1ss/wsc*ssscale
 906         v2ss=v2ss/wsc*ssscale
 907         v3ss=v3ss/wsc*ssscale
 908       else
 909         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 910       endif
 911
 912       if(me.eq.king.or..not.out1file) then
 913        write (iout,*) "Parameters of the SS-bond potential:"
 914        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
 915        " AKCT",akct
 916        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 917        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 918        write (iout,*)" HT",Ht
 919        print *,'indpdb=',indpdb,' pdbref=',pdbref
 920       endif
 921       if (indpdb.gt.0 .or. pdbref) then
 922         read(inp,'(a)') pdbfile
 923         if(me.eq.king.or..not.out1file) &
 924          write (iout,'(2a)') 'PDB data will be read from file ',&
 925          pdbfile(:ilen(pdbfile))
 926         open(ipdbin,file=pdbfile,status='old',err=33)
 927         goto 34
 928   33    write (iout,'(a)') 'Error opening PDB file.'
 929         stop
 930   34    continue
 931 !        print *,'Begin reading pdb data'
 932         call readpdb
 933 !        print *,'Finished reading pdb data'
 934         if(me.eq.king.or..not.out1file) &
 935          write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
 936          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
 937 !el        if(.not.allocated(itype_pdb))
 938         allocate(itype_pdb(nres))
 939         do i=1,nres
 940           itype_pdb(i)=itype(i,1)
 941         enddo
 942         close (ipdbin)
 943         nnt=nstart_sup
 944         nct=nstart_sup+nsup-1
 945 !el        if(.not.allocated(icont_ref))
 946         allocate(icont_ref(2,12*nres)) ! maxcont=12*maxres
 947         call contact(.false.,ncont_ref,icont_ref,co)
 948
 949         if (sideadd) then
 950          if(me.eq.king.or..not.out1file) &
 951           write(iout,*)'Adding sidechains'
 952          maxsi=1000
 953          do i=2,nres-1
 954           iti=itype(i,1)
 955           if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
 956             nsi=0
 957             fail=.true.
 958             do while (fail.and.nsi.le.maxsi)
 959               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 960               nsi=nsi+1
 961             enddo
 962             if(fail) write(iout,*)'Adding sidechain failed for res ',&
 963                     i,' after ',nsi,' trials'
 964           endif
 965          enddo
 966         endif
 967       endif
 968
 969       if (indpdb.eq.0) then
 970       nres_molec(:)=0
 971         allocate(sequence(maxres,5))
 972
 973      if (protein) then
 974 ! Read sequence if not taken from the pdb file.
 975         molec=1
 976         read (inp,*) nres_molec(molec)
 977 !        print *,'nres=',nres
 978         if (iscode.gt.0) then
 979           read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
 980         else
 981           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
 982         endif
 983 ! Convert sequence to numeric code
 984
 985         do i=1,nres_molec(molec)
 986           itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
 987         enddo
 988        endif
 989        if (nucleic) then
 990 ! Read sequence if not taken from the pdb file.
 991         molec=2
 992         read (inp,*) nres_molec(molec)
 993 !        print *,'nres=',nres
 994 !        allocate(sequence(maxres,5))
 995 !        if (iscode.gt.0) then
 996           read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
 997 ! Convert sequence to numeric code
 998
 999         do i=1,nres_molec(molec)
1000           istype(i)=sugarcode(sequence(i,molec)(1:1),i)
1001           itype(i,1)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
1002         enddo
1003        endif
1004
1005        if (ions) then
1006 ! Read sequence if not taken from the pdb file.
1007         molec=5
1008         read (inp,*) nres_molec(molec)
1009 !        print *,'nres=',nres
1010           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1011 ! Convert sequence to numeric code
1012         print *,nres_molec(molec)
1013         do i=1,nres_molec(molec)
1014           itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1015         enddo
1016        endif
1017        nres=0
1018        do i=1,5
1019         nres=nres+nres_molec(i)
1020         print *,nres_molec(i)
1021        enddo
1022
1023 ! Assign initial virtual bond lengths
1024 !elwrite(iout,*) "test_alloc"
1025         if(.not.allocated(vbld)) allocate(vbld(2*nres))
1026 !elwrite(iout,*) "test_alloc"
1027         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1028 !elwrite(iout,*) "test_alloc"
1029         do i=2,nres
1030           vbld(i)=vbl
1031           vbld_inv(i)=vblinv
1032         enddo
1033         do i=2,nres-1
1034           vbld(i+nres)=dsc(iabs(itype(i,1)))
1035           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,1)))
1036 !          write (iout,*) "i",i," itype",itype(i,1),
1037 !     &      " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1038         enddo
1039       endif
1040 !      print *,nres
1041 !      print '(20i4)',(itype(i,1),i=1,nres)
1042 !----------------------------
1043 !el reallocate tables
1044 !      do i=1,maxres2
1045 !        do j=1,3
1046 !          c_alloc(j,i)=c(j,i)
1047 !          dc_alloc(j,i)=dc(j,i)
1048 !        enddo
1049 !      enddo
1050 !      do i=1,maxres
1051 !elwrite(iout,*) "itype",i,itype(i,1)
1052 !        itype_alloc(i)=itype(i,1)
1053 !      enddo
1054
1055 !      deallocate(c)
1056 !      deallocate(dc)
1057 !      deallocate(itype)
1058 !      allocate(c(3,2*nres+4))
1059 !      allocate(dc(3,0:2*nres+2))
1060 !      allocate(itype(nres+2))
1061       allocate(itel(nres+2))
1062       itel(:)=0
1063
1064 !      do i=1,2*nres+2
1065 !        do j=1,3
1066 !          c(j,i)=c_alloc(j,i)
1067 !          dc(j,i)=dc_alloc(j,i)
1068 !        enddo
1069 !      enddo
1070 !      do i=1,nres+2
1071 !        itype(i,1)=itype_alloc(i)
1072 !        itel(i)=0
1073 !      enddo
1074 !--------------------------
1075       do i=1,nres
1076 #ifdef PROCOR
1077         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1078 #else
1079         if (itype(i,1).eq.ntyp1) then
1080 #endif
1081           itel(i)=0
1082 #ifdef PROCOR
1083         else if (iabs(itype(i+1,1)).ne.20) then
1084 #else
1085         else if (iabs(itype(i,1)).ne.20) then
1086 #endif
1087           itel(i)=1
1088         else
1089           itel(i)=2
1090         endif
1091       enddo
1092       if(me.eq.king.or..not.out1file)then
1093        write (iout,*) "ITEL"
1094        do i=1,nres-1
1095          write (iout,*) i,itype(i,1),itel(i)
1096        enddo
1097        print *,'Call Read_Bridge.'
1098       endif
1099       call read_bridge
1100 !--------------------------------
1101 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1102       call alloc_geo_arrays
1103       call alloc_ener_arrays
1104 !--------------------------------
1105 ! 8/13/98 Set limits to generating the dihedral angles
1106       do i=1,nres
1107         phibound(1,i)=-pi
1108         phibound(2,i)=pi
1109       enddo
1110       read (inp,*) ndih_constr
1111       if (ndih_constr.gt.0) then
1112         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1113         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1114
1115         read (inp,*) ftors
1116         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1117         if(me.eq.king.or..not.out1file)then
1118          write (iout,*) &
1119          'There are',ndih_constr,' constraints on phi angles.'
1120          do i=1,ndih_constr
1121           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1122          enddo
1123         endif
1124         do i=1,ndih_constr
1125           phi0(i)=deg2rad*phi0(i)
1126           drange(i)=deg2rad*drange(i)
1127         enddo
1128         if(me.eq.king.or..not.out1file) &
1129          write (iout,*) 'FTORS',ftors
1130         do i=1,ndih_constr
1131           ii = idih_constr(i)
1132           phibound(1,ii) = phi0(i)-drange(i)
1133           phibound(2,ii) = phi0(i)+drange(i)
1134         enddo
1135       endif
1136       if (with_theta_constr) then
1137 !C with_theta_constr is keyword allowing for occurance of theta constrains
1138       read (inp,*) ntheta_constr
1139 !C ntheta_constr is the number of theta constrains
1140       if (ntheta_constr.gt.0) then
1141 !C        read (inp,*) ftors
1142         allocate(itheta_constr(ntheta_constr))
1143         allocate(theta_constr0(ntheta_constr))
1144         allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1145         read (inp,*) (itheta_constr(i),theta_constr0(i), &
1146        theta_drange(i),for_thet_constr(i), &
1147        i=1,ntheta_constr)
1148 !C the above code reads from 1 to ntheta_constr
1149 !C itheta_constr(i) residue i for which is theta_constr
1150 !C theta_constr0 the global minimum value
1151 !C theta_drange is range for which there is no energy penalty
1152 !C for_thet_constr is the force constant for quartic energy penalty
1153 !C E=k*x**4
1154         if(me.eq.king.or..not.out1file)then
1155          write (iout,*) &
1156         'There are',ntheta_constr,' constraints on phi angles.'
1157          do i=1,ntheta_constr
1158           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1159          theta_drange(i), &
1160          for_thet_constr(i)
1161          enddo
1162         endif
1163         do i=1,ntheta_constr
1164           theta_constr0(i)=deg2rad*theta_constr0(i)
1165           theta_drange(i)=deg2rad*theta_drange(i)
1166         enddo
1167 !C        if(me.eq.king.or..not.out1file)
1168 !C     &   write (iout,*) 'FTORS',ftors
1169 !C        do i=1,ntheta_constr
1170 !C          ii = itheta_constr(i)
1171 !C          thetabound(1,ii) = phi0(i)-drange(i)
1172 !C          thetabound(2,ii) = phi0(i)+drange(i)
1173 !C        enddo
1174       endif ! ntheta_constr.gt.0
1175       endif! with_theta_constr
1176
1177       nnt=1
1178 #ifdef MPI
1179       if (me.eq.king) then
1180 #endif
1181        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1182        do i=1,nres
1183          write (iout,'(a3,i5,2f10.1)') &
1184          restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1185        enddo
1186 #ifdef MP
1187       endif
1188 #endif
1189       nct=nres
1190 !d      print *,'NNT=',NNT,' NCT=',NCT
1191       if (itype(1,1).eq.ntyp1) nnt=2
1192       if (itype(nres,1).eq.ntyp1) nct=nct-1
1193       if (pdbref) then
1194         if(me.eq.king.or..not.out1file) &
1195          write (iout,'(a,i3)') 'nsup=',nsup
1196         nstart_seq=nnt
1197         if (nsup.le.(nct-nnt+1)) then
1198           do i=0,nct-nnt+1-nsup
1199             if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1200               nstart_seq=nnt+i
1201               goto 111
1202             endif
1203           enddo
1204           write (iout,'(a)') &
1205                   'Error - sequences to be superposed do not match.'
1206           stop
1207         else
1208           do i=0,nsup-(nct-nnt+1)
1209             if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1210             then
1211               nstart_sup=nstart_sup+i
1212               nsup=nct-nnt+1
1213               goto 111
1214             endif
1215           enddo
1216           write (iout,'(a)') &
1217                   'Error - sequences to be superposed do not match.'
1218         endif
1219   111   continue
1220         if (nsup.eq.0) nsup=nct-nnt
1221         if (nstart_sup.eq.0) nstart_sup=nnt
1222         if (nstart_seq.eq.0) nstart_seq=nnt
1223         if(me.eq.king.or..not.out1file) &
1224          write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1225                        ' nstart_seq=',nstart_seq !,"242343453254"
1226       endif
1227 !--- Zscore rms -------
1228       if (nz_start.eq.0) nz_start=nnt
1229       if (nz_end.eq.0 .and. nsup.gt.0) then
1230         nz_end=nnt+nsup-1
1231       else if (nz_end.eq.0) then
1232         nz_end=nct
1233       endif
1234       if(me.eq.king.or..not.out1file)then
1235        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1236        write (iout,*) 'IZ_SC=',iz_sc
1237       endif
1238 !----------------------
1239       call init_int_table
1240       if (refstr) then
1241         if (.not.pdbref) then
1242           call read_angles(inp,*38)
1243           goto 39
1244    38     write (iout,'(a)') 'Error reading reference structure.'
1245 #ifdef MPI
1246           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1247           stop 'Error reading reference structure'
1248 #endif
1249    39     call chainbuild
1250           call setup_var
1251 !zscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1252           nstart_sup=nnt
1253           nstart_seq=nnt
1254           nsup=nct-nnt+1
1255           kkk=1
1256           do i=1,2*nres
1257             do j=1,3
1258               cref(j,i,kkk)=c(j,i)
1259             enddo
1260           enddo
1261           call contact(.true.,ncont_ref,icont_ref,co)
1262         endif
1263 !        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1264         call flush(iout)
1265         if (constr_dist.gt.0) call read_dist_constr
1266         write (iout,*) "After read_dist_constr nhpb",nhpb
1267         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1268         call hpb_partition
1269         if(me.eq.king.or..not.out1file) &
1270          write (iout,*) 'Contact order:',co
1271         if (pdbref) then
1272         if(me.eq.king.or..not.out1file) &
1273          write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1274         do i=1,ncont_ref
1275           do j=1,2
1276             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1277           enddo
1278           if(me.eq.king.or..not.out1file) &
1279            write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1280            icont_ref(1,i),' ',&
1281            restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1282         enddo
1283         endif
1284       endif
1285       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1286           .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1287           modecalc.ne.10) then
1288 ! If input structure hasn't been supplied from the PDB file read or generate
1289 ! initial geometry.
1290         if (iranconf.eq.0 .and. .not. extconf) then
1291           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1292            write (iout,'(a)') 'Initial geometry will be read in.'
1293           if (read_cart) then
1294             read(inp,'(8f10.5)',end=36,err=36) &
1295              ((c(l,k),l=1,3),k=1,nres),&
1296              ((c(l,k+nres),l=1,3),k=nnt,nct)
1297             write (iout,*) "Exit READ_CART"
1298             write (iout,'(8f10.5)') &
1299              ((c(l,k),l=1,3),k=1,nres),&
1300              ((c(l,k+nres),l=1,3),k=nnt,nct)
1301             call int_from_cart1(.true.)
1302             write (iout,*) "Finish INT_TO_CART"
1303             do i=1,nres-1
1304               do j=1,3
1305                 dc(j,i)=c(j,i+1)-c(j,i)
1306                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1307               enddo
1308             enddo
1309             do i=nnt,nct
1310               if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1311                 do j=1,3
1312                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1313                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1314                 enddo
1315               endif
1316             enddo
1317             return
1318           else
1319             call read_angles(inp,*36)
1320           endif
1321           goto 37
1322    36     write (iout,'(a)') 'Error reading angle file.'
1323 #ifdef MPI
1324           call mpi_finalize( MPI_COMM_WORLD,IERR )
1325 #endif
1326           stop 'Error reading angle file.'
1327    37     continue
1328         else if (extconf) then
1329          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1330           write (iout,'(a)') 'Extended chain initial geometry.'
1331          do i=3,nres
1332           theta(i)=90d0*deg2rad
1333          enddo
1334          do i=4,nres
1335           phi(i)=180d0*deg2rad
1336          enddo
1337          do i=2,nres-1
1338           alph(i)=110d0*deg2rad
1339          enddo
1340 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1341          do i=2,nres-1
1342           omeg(i)=-120d0*deg2rad
1343           if (itype(i,1).le.0) omeg(i)=-omeg(i)
1344          enddo
1345         else
1346           if(me.eq.king.or..not.out1file) &
1347            write (iout,'(a)') 'Random-generated initial geometry.'
1348
1349
1350 #ifdef MPI
1351           if (me.eq.king  .or. fg_rank.eq.0 .and. &
1352                  ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1353 #endif
1354             do itrial=1,100
1355               itmp=1
1356               call gen_rand_conf(itmp,*30)
1357               goto 40
1358    30         write (iout,*) 'Failed to generate random conformation',&
1359                 ', itrial=',itrial
1360               write (*,*) 'Processor:',me,&
1361                 ' Failed to generate random conformation',&
1362                 ' itrial=',itrial
1363               call intout
1364
1365 #ifdef AIX
1366               call flush_(iout)
1367 #else
1368               call flush(iout)
1369 #endif
1370             enddo
1371             write (iout,'(a,i3,a)') 'Processor:',me,&
1372               ' error in generating random conformation.'
1373             write (*,'(a,i3,a)') 'Processor:',me,&
1374               ' error in generating random conformation.'
1375             call flush(iout)
1376 #ifdef MPI
1377             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1378    40       continue
1379           endif
1380 #else
1381           do itrial=1,100
1382             itmp=1
1383             call gen_rand_conf(itmp,*335)
1384             goto 40
1385   335       write (iout,*) 'Failed to generate random conformation',&
1386               ', itrial=',itrial
1387             write (*,*) 'Failed to generate random conformation',&
1388               ', itrial=',itrial
1389           enddo
1390           write (iout,'(a,i3,a)') 'Processor:',me,&
1391             ' error in generating random conformation.'
1392           write (*,'(a,i3,a)') 'Processor:',me,&
1393             ' error in generating random conformation.'
1394           stop
1395    40     continue
1396 #endif
1397         endif
1398       elseif (modecalc.eq.4) then
1399         read (inp,'(a)') intinname
1400         open (intin,file=intinname,status='old',err=333)
1401         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1402         write (iout,'(a)') 'intinname',intinname
1403         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1404         goto 334
1405   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1406 #ifdef MPI
1407         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1408 #endif
1409         stop 'Error opening angle file.'
1410   334   continue
1411
1412       endif
1413 ! Generate distance constraints, if the PDB structure is to be regularized.
1414       if (nthread.gt.0) then
1415         call read_threadbase
1416       endif
1417       call setup_var
1418 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1419       if (me.eq.king .or. .not. out1file) &
1420        call intout
1421 !elwrite (iout,*)"alph(i)*rad2deg",(alph(i)*rad2deg, i=1,nres)
1422       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1423         write (iout,'(/a,i3,a)') &
1424         'The chain contains',ns,' disulfide-bridging cysteines.'
1425         write (iout,'(20i4)') (iss(i),i=1,ns)
1426        if (dyn_ss) then
1427           write(iout,*)"Running with dynamic disulfide-bond formation"
1428        else
1429         write (iout,'(/a/)') 'Pre-formed links are:'
1430         do i=1,nss
1431           i1=ihpb(i)-nres
1432           i2=jhpb(i)-nres
1433           it1=itype(i1,1)
1434           it2=itype(i2,1)
1435           if (me.eq.king.or..not.out1file) &
1436           write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1437           restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1438           ebr,forcon(i)
1439         enddo
1440         write (iout,'(a)')
1441        endif
1442       endif
1443       if (ns.gt.0.and.dyn_ss) then
1444           do i=nss+1,nhpb
1445             ihpb(i-nss)=ihpb(i)
1446             jhpb(i-nss)=jhpb(i)
1447             forcon(i-nss)=forcon(i)
1448             dhpb(i-nss)=dhpb(i)
1449           enddo
1450           nhpb=nhpb-nss
1451           nss=0
1452           call hpb_partition
1453           do i=1,ns
1454             dyn_ss_mask(iss(i))=.true.
1455           enddo
1456       endif
1457       if (i2ndstr.gt.0) call secstrp2dihc
1458 !      call geom_to_var(nvar,x)
1459 !      call etotal(energia(0))
1460 !      call enerprint(energia(0))
1461 !      call briefout(0,etot)
1462 !      stop
1463 !d    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1464 !d    write (iout,'(a)') 'Variable list:'
1465 !d    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1466 #ifdef MPI
1467       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1468         write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1469         'Processor',myrank,': end reading molecular data.'
1470 #endif
1471       return
1472       end subroutine molread
1473 !-----------------------------------------------------------------------------
1474 !-----------------------------------------------------------------------------
1475       end module io