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