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,molnum(i))
 229                  ii = ii+1
 230                enddo
 231              endif
 232            enddo
 233            call chainbuild
 234            call etotal(energia)
 235            fail = (energia(0).ge.1.0d20)
 236            icount_fail=icount_fail+1
 237
 238            ENDDO
 239
 240            if (icount_fail.gt.maxcount_fail) then
 241              write (iout,*) &
 242              'Failed to generate non-overlaping near-native conf.',&
 243              m
 244            endif
 245
 246            do j=2,nres-1
 247              dihang_in(1,j,1,m)=theta(j+1)
 248              dihang_in(2,j,1,m)=phi(j+2)
 249              dihang_in(3,j,1,m)=alph(j)
 250              dihang_in(4,j,1,m)=omeg(j)
 251            enddo
 252            dihang_in(2,nres-1,1,m)=0.0d0
 253          enddo
 254
 255 !      do m=1,n
 256 !        write(icsa_native_int,*) m,e
 257 !         write(icsa_native_int,200) (dihang_in(1,k,1,m),k=2,nres-1)
 258 !         write(icsa_native_int,200) (dihang_in(2,k,1,m),k=2,nres-2)
 259 !         write(icsa_native_int,200) (dihang_in(3,k,1,m),k=2,nres-1)
 260 !         write(icsa_native_int,200) (dihang_in(4,k,1,m),k=2,nres-1)
 261 !      enddo
 262 !     close(icsa_native_int)
 263
 264 !      do m=mm+2,n
 265 !       do i=1,4
 266 !        do j=2,nres-1
 267 !         dihang_in(i,j,1,m)=dihang_in(i,j,1,m)*deg2rad
 268 !        enddo
 269 !       enddo
 270 !      enddo
 271
 272       call dihang_to_c(dihang_in(1,1,1,1))
 273
 274 ! Store c to cref (they are in COMMON.CHAIN).
 275       do k=1,2*nres
 276        do kk=1,3
 277         crefjlee(kk,k)=c(kk,k)
 278        enddo
 279       enddo
 280
 281       call contact(.true.,ncont_ref,icont_ref,co)
 282
 283 !      do k=1,nres
 284 !       write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
 285 !      enddo
 286       close(icsa_native_int)
 287
 288   200 format (8f10.4)
 289
 290       return
 291       end subroutine from_int
 292 !-----------------------------------------------------------------------------
 293       subroutine dihang_to_c(aarray)
 294
 295       use geometry_data
 296       use csa_data
 297       use geometry, only:chainbuild
 298 !      implicit real*8 (a-h,o-z)
 299 !      include 'DIMENSIONS'
 300 !      include 'COMMON.CSA'
 301 !      include 'COMMON.BANK'
 302 !      include 'COMMON.CHAIN'
 303 !      include 'COMMON.GEO'
 304 !      include 'COMMON.VAR'
 305       integer :: i
 306       real(kind=8),dimension(mxang,nres,mxch) :: aarray !(mxang,maxres,mxch)
 307
 308 !     do i=4,nres
 309 !      phi(i)=dihang_in(1,i-2,1,1)
 310 !     enddo
 311       do i=2,nres-1
 312        theta(i+1)=aarray(1,i,1)
 313        phi(i+2)=aarray(2,i,1)
 314        alph(i)=aarray(3,i,1)
 315        omeg(i)=aarray(4,i,1)
 316       enddo
 317
 318       call chainbuild
 319
 320       return
 321       end subroutine dihang_to_c
 322 !-----------------------------------------------------------------------------
 323 ! geomout.F
 324 !-----------------------------------------------------------------------------
 325 #ifdef NOXDR
 326       subroutine cartout(time)
 327 #else
 328       subroutine cartoutx(time)
 329 #endif
 330       use geometry_data, only: c,nres
 331       use energy_data
 332       use MD_data, only: potE,t_bath
 333 !      implicit real*8 (a-h,o-z)
 334 !      include 'DIMENSIONS'
 335 !      include 'COMMON.CHAIN'
 336 !      include 'COMMON.INTERACT'
 337 !      include 'COMMON.NAMES'
 338 !      include 'COMMON.IOUNITS'
 339 !      include 'COMMON.HEADER'
 340 !      include 'COMMON.SBRIDGE'
 341 !      include 'COMMON.DISTFIT'
 342 !      include 'COMMON.MD'
 343       real(kind=8) :: time
 344 !el  local variables
 345       integer :: j,k,i
 346
 347 #if defined(AIX) || defined(PGI)
 348       open(icart,file=cartname,position="append")
 349 #else
 350       open(icart,file=cartname,access="append")
 351 #endif
 352       write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 353       if (dyn_ss) then
 354        write (icart,'(i4,$)') &
 355          nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)
 356       else
 357        write (icart,'(i4,$)') &
 358          nss,(ihpb(j),jhpb(j),j=1,nss)
 359        endif
 360        write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,&
 361        (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),&
 362        (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 363       write (icart,'(8f10.5)') &
 364        ((c(k,j),k=1,3),j=1,nres),&
 365        ((c(k,j+nres),k=1,3),j=nnt,nct)
 366       close(icart)
 367       return
 368
 369 #ifdef NOXDR
 370       end subroutine cartout
 371 #else
 372       end subroutine cartoutx
 373 #endif
 374 !-----------------------------------------------------------------------------
 375 #ifndef NOXDR
 376       subroutine cartout(time)
 377 !      implicit real*8 (a-h,o-z)
 378 !      include 'DIMENSIONS'
 379       use geometry_data, only: c,nres
 380       use energy_data
 381       use MD_data, only: potE,t_bath
 382 #ifdef MPI
 383       use MPI_data
 384       include 'mpif.h'
 385 !      include 'COMMON.SETUP'
 386 #else
 387       integer,parameter :: me=0
 388 #endif
 389 !      include 'COMMON.CHAIN'
 390 !      include 'COMMON.INTERACT'
 391 !      include 'COMMON.NAMES'
 392 !      include 'COMMON.IOUNITS'
 393 !      include 'COMMON.HEADER'
 394 !      include 'COMMON.SBRIDGE'
 395 !      include 'COMMON.DISTFIT'
 396 !      include 'COMMON.MD'
 397       real(kind=8) :: time
 398       integer :: iret,itmp
 399       real(kind=4) :: prec
 400       real(kind=4),dimension(3,2*nres+2) :: xcoord      !(3,maxres2+2)  (maxres2=2*maxres
 401 !el  local variables
 402       integer :: j,i,ixdrf
 403
 404 #ifdef AIX
 405       call xdrfopen_(ixdrf,cartname, "a", iret)
 406       call xdrffloat_(ixdrf, real(time), iret)
 407       call xdrffloat_(ixdrf, real(potE), iret)
 408       call xdrffloat_(ixdrf, real(uconst), iret)
 409       call xdrffloat_(ixdrf, real(uconst_back), iret)
 410       call xdrffloat_(ixdrf, real(t_bath), iret)
 411       call xdrfint_(ixdrf, nss, iret)
 412       do j=1,nss
 413        if (dyn_ss) then
 414         call xdrfint_(ixdrf, idssb(j)+nres, iret)
 415         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 416        else
 417         call xdrfint_(ixdrf, ihpb(j), iret)
 418         call xdrfint_(ixdrf, jhpb(j), iret)
 419        endif
 420       enddo
 421       call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
 422       do i=1,nfrag
 423         call xdrffloat_(ixdrf, real(qfrag(i)), iret)
 424       enddo
 425       do i=1,npair
 426         call xdrffloat_(ixdrf, real(qpair(i)), iret)
 427       enddo
 428       do i=1,nfrag_back
 429         call xdrffloat_(ixdrf, real(utheta(i)), iret)
 430         call xdrffloat_(ixdrf, real(ugamma(i)), iret)
 431         call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
 432       enddo
 433 #else
 434       call xdrfopen(ixdrf,cartname, "a", iret)
 435       call xdrffloat(ixdrf, real(time), iret)
 436       call xdrffloat(ixdrf, real(potE), iret)
 437       call xdrffloat(ixdrf, real(uconst), iret)
 438       call xdrffloat(ixdrf, real(uconst_back), iret)
 439       call xdrffloat(ixdrf, real(t_bath), iret)
 440       call xdrfint(ixdrf, nss, iret)
 441       do j=1,nss
 442        if (dyn_ss) then
 443         call xdrfint(ixdrf, idssb(j)+nres, iret)
 444         call xdrfint(ixdrf, jdssb(j)+nres, iret)
 445        else
 446         call xdrfint(ixdrf, ihpb(j), iret)
 447         call xdrfint(ixdrf, jhpb(j), iret)
 448        endif
 449       enddo
 450       call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
 451       do i=1,nfrag
 452         call xdrffloat(ixdrf, real(qfrag(i)), iret)
 453       enddo
 454       do i=1,npair
 455         call xdrffloat(ixdrf, real(qpair(i)), iret)
 456       enddo
 457       do i=1,nfrag_back
 458         call xdrffloat(ixdrf, real(utheta(i)), iret)
 459         call xdrffloat(ixdrf, real(ugamma(i)), iret)
 460         call xdrffloat(ixdrf, real(uscdiff(i)), iret)
 461       enddo
 462 #endif
 463       prec=10000.0
 464       do i=1,nres
 465        do j=1,3
 466         xcoord(j,i)=c(j,i)
 467        enddo
 468       enddo
 469       do i=nnt,nct
 470        do j=1,3
 471         xcoord(j,nres+i-nnt+1)=c(j,i+nres)
 472        enddo
 473       enddo
 474
 475       itmp=nres+nct-nnt+1
 476 #ifdef AIX
 477       call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
 478       call xdrfclose_(ixdrf, iret)
 479 #else
 480       call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
 481       call xdrfclose(ixdrf, iret)
 482 #endif
 483       return
 484       end subroutine cartout
 485 #endif
 486 !-----------------------------------------------------------------------------
 487       subroutine statout(itime)
 488
 489       use energy_data
 490       use control_data
 491       use MD_data
 492       use MPI_data
 493       use compare, only:rms_nac_nnc
 494 !      implicit real*8 (a-h,o-z)
 495 !      include 'DIMENSIONS'
 496 !      include 'COMMON.CONTROL'
 497 !      include 'COMMON.CHAIN'
 498 !      include 'COMMON.INTERACT'
 499 !      include 'COMMON.NAMES'
 500 !      include 'COMMON.IOUNITS'
 501 !      include 'COMMON.HEADER'
 502 !      include 'COMMON.SBRIDGE'
 503 !      include 'COMMON.DISTFIT'
 504 !      include 'COMMON.MD'
 505 !      include 'COMMON.REMD'
 506 !      include 'COMMON.SETUP'
 507       integer :: itime
 508       real(kind=8),dimension(0:n_ene) :: energia
 509 !      double precision gyrate
 510 !el      external gyrate
 511 !el      common /gucio/ cm
 512       character(len=256) :: line1,line2
 513       character(len=4) :: format1,format2
 514       character(len=30) :: format
 515 !el  local variables
 516       integer :: i,ii1,ii2,j
 517       real(kind=8) :: rms,frac,frac_nn,co,distance
 518
 519 #ifdef AIX
 520       if(itime.eq.0) then
 521        open(istat,file=statname,position="append")
 522       endif
 523 #else
 524 #ifdef PGI
 525       open(istat,file=statname,position="append")
 526 #else
 527       open(istat,file=statname,access="append")
 528 #endif
 529 #endif
 530        if (AFMlog.gt.0) then
 531        if (refstr) then
 532          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 533           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')&
 534                itime,totT,EK,potE,totE,&
 535                rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
 536                potEcomp(23),me
 537           format1="a133"
 538          else
 539 !C          print *,'A CHUJ',potEcomp(23)
 540           write (line1,'(i10,f15.2,7f12.3,i5,$)') &
 541                 itime,totT,EK,potE,totE,&
 542                 kinetic_T,t_bath,gyrate(),&
 543                 potEcomp(23),me
 544           format1="a114"
 545         endif
 546        else if (selfguide.gt.0) then
 547        distance=0.0
 548        do j=1,3
 549        distance=distance+(c(j,afmend)-c(j,afmbeg))**2
 550        enddo
 551        distance=dsqrt(distance)
 552        if (refstr) then
 553          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 554           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
 555          f9.3,i5,$)') &
 556                itime,totT,EK,potE,totE,&
 557                rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
 558                distance,potEcomp(23),me
 559           format1="a133"
 560 !C          print *,"CHUJOWO"
 561          else
 562 !C          print *,'A CHUJ',potEcomp(23)
 563           write (line1,'(i10,f15.2,8f12.3,i5,$)')&
 564                 itime,totT,EK,potE,totE, &
 565                 kinetic_T,t_bath,gyrate(),&
 566                 distance,potEcomp(23),me
 567           format1="a114"
 568         endif
 569        else
 570        if (refstr) then
 571          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 572           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
 573                 itime,totT,EK,potE,totE,&
 574                 rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
 575           format1="a133"
 576         else
 577           write (line1,'(i10,f15.2,7f12.3,i5,$)') &
 578                  itime,totT,EK,potE,totE,&
 579                  amax,kinetic_T,t_bath,gyrate(),me
 580           format1="a114"
 581         endif
 582         ENDIF
 583         if(usampl.and.totT.gt.eq_time) then
 584            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,&
 585             (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),&
 586             (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 587            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair &
 588                    +21*nfrag_back
 589         else
 590            format2="a001"
 591            line2=' '
 592         endif
 593         if (print_compon) then
 594           if(itime.eq.0) then
 595            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 596                                                            ",20a12)"
 597            write (istat,format) "#","",&
 598             (ename(print_order(i)),i=1,nprint_ene)
 599           endif
 600           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 601                                                            ",20f12.3)"
 602           write (istat,format) line1,line2,&
 603             (potEcomp(print_order(i)),i=1,nprint_ene)
 604         else
 605           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
 606           write (istat,format) line1,line2
 607         endif
 608 #if defined(AIX)
 609         call flush(istat)
 610 #else
 611         close(istat)
 612 #endif
 613       return
 614       end subroutine  statout
 615 !-----------------------------------------------------------------------------
 616 ! readrtns_CSA.F
 617 !-----------------------------------------------------------------------------
 618       subroutine readrtns
 619
 620       use control_data
 621       use energy_data
 622       use MPI_data
 623       use muca_md, only:read_muca
 624 !      implicit real*8 (a-h,o-z)
 625 !      include 'DIMENSIONS'
 626 #ifdef MPI
 627       include 'mpif.h'
 628 #endif
 629 !      include 'COMMON.SETUP'
 630 !      include 'COMMON.CONTROL'
 631 !      include 'COMMON.SBRIDGE'
 632 !      include 'COMMON.IOUNITS'
 633       logical :: file_exist
 634       integer :: i
 635 ! Read force-field parameters except weights
 636 !      call parmread
 637 ! Read job setup parameters
 638       call read_control
 639 ! Read force-field parameters except weights
 640       call parmread
 641
 642 ! Read control parameters for energy minimzation if required
 643       if (minim) call read_minim
 644 ! Read MCM control parameters if required
 645       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
 646 ! Read MD control parameters if reqjuired
 647       if (modecalc.eq.12) call read_MDpar
 648 ! Read MREMD control parameters if required
 649       if (modecalc.eq.14) then
 650          call read_MDpar
 651          call read_REMDpar
 652       endif
 653 ! Read MUCA control parameters if required
 654       if (lmuca) call read_muca
 655 ! Read CSA control parameters if required (from fort.40 if exists
 656 ! otherwise from general input file)
 657       if (modecalc.eq.8) then
 658        inquire (file="fort.40",exist=file_exist)
 659        if (.not.file_exist) call csaread
 660       endif
 661 !fmc      if (modecalc.eq.10) call mcmfread
 662 ! Read molecule information, molecule geometry, energy-term weights, and
 663 ! restraints if requested
 664       call molread
 665 ! Print restraint information
 666 #ifdef MPI
 667       if (.not. out1file .or. me.eq.king) then
 668 #endif
 669       if (nhpb.gt.nss) &
 670       write (iout,'(a,i5,a)') "The following",nhpb-nss,&
 671        " distance constraints have been imposed"
 672       do i=nss+1,nhpb
 673         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
 674       enddo
 675 #ifdef MPI
 676       endif
 677 #endif
 678 !      print *,"Processor",myrank," leaves READRTNS"
 679 !      write(iout,*) "end readrtns"
 680       return
 681       end subroutine readrtns
 682 !-----------------------------------------------------------------------------
 683       subroutine molread
 684 !
 685 ! Read molecular data.
 686 !
 687 !      use control, only: ilen
 688       use control_data
 689       use geometry_data
 690       use energy_data
 691       use energy
 692       use compare_data
 693       use MD_data, only: t_bath
 694       use MPI_data
 695       use compare, only:seq_comp,contact
 696       use control
 697 !      implicit real*8 (a-h,o-z)
 698 !      include 'DIMENSIONS'
 699 #ifdef MPI
 700       include 'mpif.h'
 701       integer :: error_msg,ierror,ierr,ierrcode
 702 #endif
 703 !      include 'COMMON.IOUNITS'
 704 !      include 'COMMON.GEO'
 705 !      include 'COMMON.VAR'
 706 !      include 'COMMON.INTERACT'
 707 !      include 'COMMON.LOCAL'
 708 !      include 'COMMON.NAMES'
 709 !      include 'COMMON.CHAIN'
 710 !      include 'COMMON.FFIELD'
 711 !      include 'COMMON.SBRIDGE'
 712 !      include 'COMMON.HEADER'
 713 !      include 'COMMON.CONTROL'
 714 !      include 'COMMON.DBASE'
 715 !      include 'COMMON.THREAD'
 716 !      include 'COMMON.CONTACTS'
 717 !      include 'COMMON.TORCNSTR'
 718 !      include 'COMMON.TIME1'
 719 !      include 'COMMON.BOUNDS'
 720 !      include 'COMMON.MD'
 721 !      include 'COMMON.SETUP'
 722       character(len=4),dimension(:,:),allocatable :: sequence   !(maxres,maxmolecules)
 723 !      integer :: rescode
 724 !      double precision x(maxvar)
 725       character(len=256) :: pdbfile
 726       character(len=800) :: weightcard
 727       character(len=80) :: weightcard_t!,ucase
 728 !      integer,dimension(:),allocatable :: itype_pdb    !(maxres)
 729 !      common /pizda/ itype_pdb
 730       logical :: fail   !seq_comp,
 731       real(kind=8) :: energia(0:n_ene)
 732 !      integer ilen
 733 !el      external ilen
 734 !el local varables
 735       integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
 736
 737       real(kind=8),dimension(3,maxres2+2) :: c_alloc
 738       real(kind=8),dimension(3,0:maxres2) :: dc_alloc
 739       real(kind=8),dimension(:,:), allocatable :: secprob
 740       integer,dimension(maxres) :: itype_alloc
 741
 742       integer :: iti,nsi,maxsi,itrial,itmp
 743       real(kind=8) :: wlong,scalscp,co,ssscale,phihel,phibet,sigmahel,&
 744       sigmabet,sumv,nres_temp
 745       allocate(weights(n_ene))
 746 !-----------------------------
 747       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
 748       allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
 749       allocate(itype(maxres,5)) !(maxres)
 750       allocate(istype(maxres))
 751 !
 752 ! Zero out tables.
 753 !
 754       c(:,:)=0.0D0
 755       dc(:,:)=0.0D0
 756       itype(:,:)=0
 757 !-----------------------------
 758 !
 759 ! Body
 760 !
 761 ! Read weights of the subsequent energy terms.
 762       call card_concat(weightcard,.true.)
 763       call reada(weightcard,'WLONG',wlong,1.0D0)
 764       call reada(weightcard,'WSC',wsc,wlong)
 765       call reada(weightcard,'WSCP',wscp,wlong)
 766       call reada(weightcard,'WELEC',welec,1.0D0)
 767       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 768       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 769       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 770       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 771       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 772       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 773       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 774       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 775       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 776       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 777       call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
 778       call reada(weightcard,'WELPP',welpp,0.0d0)
 779       call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
 780       call reada(weightcard,'WELPSB',welpsb,0.0D0)
 781       call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
 782       call reada(weightcard,'WELSB',welsb,0.0D0)
 783       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
 784       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
 785       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
 786       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
 787 !      print *,"KUR..",wtor_nucl
 788       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
 789       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
 790       call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
 791       call reada(weightcard,'WBOND',wbond,1.0D0)
 792       call reada(weightcard,'WTOR',wtor,1.0D0)
 793       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 794       call reada(weightcard,'WSHIELD',wshield,0.05D0)
 795       call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
 796       call reada(weightcard,'WLT',wliptran,1.0D0)
 797       call reada(weightcard,'WTUBE',wtube,1.0d0)
 798       call reada(weightcard,'WANG',wang,1.0D0)
 799       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 800       call reada(weightcard,'SCAL14',scal14,0.4D0)
 801       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 802       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 803       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 804       call reada(weightcard,'TEMP0',temp0,300.0d0)
 805       call reada(weightcard,'WSCBASE',wscbase,0.0D0)
 806       if (index(weightcard,'SOFT').gt.0) ipot=6
 807       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
 808       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
 809       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
 810       call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
 811       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
 812       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
 813       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
 814
 815 ! 12/1/95 Added weight for the multi-body term WCORR
 816       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 817       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 818       weights(1)=wsc
 819       weights(2)=wscp
 820       weights(3)=welec
 821       weights(4)=wcorr
 822       weights(5)=wcorr5
 823       weights(6)=wcorr6
 824       weights(7)=wel_loc
 825       weights(8)=wturn3
 826       weights(9)=wturn4
 827       weights(10)=wturn6
 828       weights(11)=wang
 829       weights(12)=wscloc
 830       weights(13)=wtor
 831       weights(14)=wtor_d
 832       weights(15)=wstrain
 833       weights(16)=wvdwpp
 834       weights(17)=wbond
 835       weights(18)=scal14
 836       weights(21)=wsccor
 837           weights(26)=wvdwpp_nucl
 838           weights(27)=welpp
 839           weights(28)=wvdwpsb
 840           weights(29)=welpsb
 841           weights(30)=wvdwsb
 842           weights(31)=welsb
 843           weights(32)=wbond_nucl
 844           weights(33)=wang_nucl
 845           weights(34)=wsbloc
 846           weights(35)=wtor_nucl
 847           weights(36)=wtor_d_nucl
 848           weights(37)=wcorr_nucl
 849           weights(38)=wcorr3_nucl
 850           weights(41)=wcatcat
 851           weights(42)=wcatprot
 852           weights(46)=wscbase
 853           weights(47)=wpepbase
 854           weights(48)=wscpho
 855           weights(49)=wpeppho
 856           weights(50)=wcatnucl
 857
 858       if(me.eq.king.or..not.out1file) &
 859        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 860         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 861         wturn4,wturn6
 862    10 format (/'Energy-term weights (unscaled):'// &
 863        'WSCC=   ',f10.6,' (SC-SC)'/ &
 864        'WSCP=   ',f10.6,' (SC-p)'/ &
 865        'WELEC=  ',f10.6,' (p-p electr)'/ &
 866        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 867        'WBOND=  ',f10.6,' (stretching)'/ &
 868        'WANG=   ',f10.6,' (bending)'/ &
 869        'WSCLOC= ',f10.6,' (SC local)'/ &
 870        'WTOR=   ',f10.6,' (torsional)'/ &
 871        'WTORD=  ',f10.6,' (double torsional)'/ &
 872        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 873        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 874        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 875        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 876        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 877        'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
 878        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 879        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 880        'WTURN6= ',f10.6,' (turns, 6th order)')
 881       if(me.eq.king.or..not.out1file)then
 882        if (wcorr4.gt.0.0d0) then
 883         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
 884          'between contact pairs of peptide groups'
 885         write (iout,'(2(a,f5.3/))') &
 886         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
 887         'Range of quenching the correlation terms:',2*delt_corr
 888        else if (wcorr.gt.0.0d0) then
 889         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
 890          'between contact pairs of peptide groups'
 891        endif
 892        write (iout,'(a,f8.3)') &
 893         'Scaling factor of 1,4 SC-p interactions:',scal14
 894        write (iout,'(a,f8.3)') &
 895         'General scaling factor of SC-p interactions:',scalscp
 896       endif
 897       r0_corr=cutoff_corr-delt_corr
 898       do i=1,ntyp
 899         aad(i,1)=scalscp*aad(i,1)
 900         aad(i,2)=scalscp*aad(i,2)
 901         bad(i,1)=scalscp*bad(i,1)
 902         bad(i,2)=scalscp*bad(i,2)
 903       enddo
 904       call rescale_weights(t_bath)
 905       if(me.eq.king.or..not.out1file) &
 906        write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 907         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 908         wturn4,wturn6
 909    22 format (/'Energy-term weights (scaled):'// &
 910        'WSCC=   ',f10.6,' (SC-SC)'/ &
 911        'WSCP=   ',f10.6,' (SC-p)'/ &
 912        'WELEC=  ',f10.6,' (p-p electr)'/ &
 913        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 914        'WBOND=  ',f10.6,' (stretching)'/ &
 915        'WANG=   ',f10.6,' (bending)'/ &
 916        'WSCLOC= ',f10.6,' (SC local)'/ &
 917        'WTOR=   ',f10.6,' (torsional)'/ &
 918        'WTORD=  ',f10.6,' (double torsional)'/ &
 919        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 920        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 921        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 922        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 923        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 924        'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
 925        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 926        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 927        'WTURN6= ',f10.6,' (turns, 6th order)')
 928       if(me.eq.king.or..not.out1file) &
 929        write (iout,*) "Reference temperature for weights calculation:",&
 930         temp0
 931       call reada(weightcard,"D0CM",d0cm,3.78d0)
 932       call reada(weightcard,"AKCM",akcm,15.1d0)
 933       call reada(weightcard,"AKTH",akth,11.0d0)
 934       call reada(weightcard,"AKCT",akct,12.0d0)
 935       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 936       call reada(weightcard,"V2SS",v2ss,7.61d0)
 937       call reada(weightcard,"V3SS",v3ss,13.7d0)
 938       call reada(weightcard,"EBR",ebr,-5.50D0)
 939       call reada(weightcard,"ATRISS",atriss,0.301D0)
 940       call reada(weightcard,"BTRISS",btriss,0.021D0)
 941       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 942       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 943       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
 944       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 945
 946       call reada(weightcard,"HT",Ht,0.0D0)
 947       if (dyn_ss) then
 948        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
 949         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
 950         akcm=akcm/wsc*ssscale
 951         akth=akth/wsc*ssscale
 952         akct=akct/wsc*ssscale
 953         v1ss=v1ss/wsc*ssscale
 954         v2ss=v2ss/wsc*ssscale
 955         v3ss=v3ss/wsc*ssscale
 956       else
 957         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 958       endif
 959
 960       if(me.eq.king.or..not.out1file) then
 961        write (iout,*) "Parameters of the SS-bond potential:"
 962        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
 963        " AKCT",akct
 964        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 965        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 966        write (iout,*)" HT",Ht
 967        print *,'indpdb=',indpdb,' pdbref=',pdbref
 968       endif
 969       if (indpdb.gt.0 .or. pdbref) then
 970         read(inp,'(a)') pdbfile
 971         if(me.eq.king.or..not.out1file) &
 972          write (iout,'(2a)') 'PDB data will be read from file ',&
 973          pdbfile(:ilen(pdbfile))
 974         open(ipdbin,file=pdbfile,status='old',err=33)
 975         goto 34
 976   33    write (iout,'(a)') 'Error opening PDB file.'
 977         stop
 978   34    continue
 979 !        print *,'Begin reading pdb data'
 980         call readpdb
 981         if (.not.allocated(crefjlee)) allocate (crefjlee(3,2*nres+2))
 982         do i=1,2*nres
 983           do j=1,3
 984             crefjlee(j,i)=c(j,i)
 985           enddo
 986         enddo
 987 #ifdef DEBUG
 988         do i=1,nres
 989           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 990      &      (crefjlee(j,i+nres),j=1,3)
 991         enddo
 992 #endif
 993
 994 !        call int_from_cart1(.true.)
 995
 996 !        print *,'Finished reading pdb data'
 997         if(me.eq.king.or..not.out1file) &
 998          write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
 999          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
1000 !el        if(.not.allocated(itype_pdb))
1001         allocate(itype_pdb(nres))
1002         do i=1,nres
1003           itype_pdb(i)=itype(i,1)
1004         enddo
1005         close (ipdbin)
1006         nnt=nstart_sup
1007         nct=nstart_sup+nsup-1
1008 !el        if(.not.allocated(icont_ref))
1009         allocate(icont_ref(2,(nres/2)*nres)) ! maxcont=12*maxres
1010         call contact(.false.,ncont_ref,icont_ref,co)
1011
1012         if (sideadd) then
1013          if(me.eq.king.or..not.out1file) &
1014           write(iout,*)'Adding sidechains'
1015          maxsi=1000
1016          do i=2,nres-1
1017           iti=itype(i,1)
1018           if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
1019             nsi=0
1020             fail=.true.
1021             do while (fail.and.nsi.le.maxsi)
1022               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i))
1023               nsi=nsi+1
1024             enddo
1025             if(fail) write(iout,*)'Adding sidechain failed for res ',&
1026                     i,' after ',nsi,' trials'
1027           endif
1028          enddo
1029         endif
1030       endif
1031       print *,"CZY TU DOCHODZE"
1032       if (indpdb.eq.0) then
1033       nres_molec(:)=0
1034         allocate(sequence(maxres,5))
1035 !      itype(:,:)=0
1036       itmp=0
1037       if (protein) then
1038 ! Read sequence if not taken from the pdb file.
1039         molec=1
1040         read (inp,*) nres_molec(molec)
1041         print *,'nres=',nres
1042         if (iscode.gt.0) then
1043           read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
1044         else
1045           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1046         endif
1047 !        read(inp,*) weightcard_t
1048 !        print *,"po seq" weightcard_t
1049 ! Convert sequence to numeric code
1050
1051         do i=1,nres_molec(molec)
1052           itmp=itmp+1
1053           itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1054           print *,itype(i,1)
1055
1056         enddo
1057        endif
1058 !        read(inp,*) weightcard_t
1059 !        print *,"po seq", weightcard_t
1060
1061        if (nucleic) then
1062 ! Read sequence if not taken from the pdb file.
1063         molec=2
1064         read (inp,*) nres_molec(molec)
1065         print *,'nres=',nres_molec(molec)
1066 !        allocate(sequence(maxres,5))
1067 !        if (iscode.gt.0) then
1068           read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
1069           print *,"KUR**"
1070           print *,(sequence(i,molec),i=1,nres_molec(molec))
1071 ! Convert sequence to numeric code
1072
1073         do i=1,nres_molec(molec)
1074           itmp=itmp+1
1075           istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)
1076           sequence(i,molec)=sequence(i,molec)(1:2)
1077           itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1078           write(iout,*) "NUCLE=", itype(itmp,molec)
1079         enddo
1080        endif
1081
1082        if (ions) then
1083 ! Read sequence if not taken from the pdb file.
1084         molec=5
1085         read (inp,*) nres_molec(molec)
1086 !        print *,'nres=',nres
1087           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1088 ! Convert sequence to numeric code
1089         print *,nres_molec(molec)
1090         do i=1,nres_molec(molec)
1091           itmp=itmp+1
1092           print *,itmp,"itmp"
1093           itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1094         enddo
1095        endif
1096        nres=0
1097        do i=1,5
1098         nres=nres+nres_molec(i)
1099         print *,"nres_molec",nres,nres_molec(i)
1100        enddo
1101
1102 ! Assign initial virtual bond lengths
1103         if(.not.allocated(molnum)) then
1104          allocate(molnum(nres+1))
1105          itmp=0
1106         do i=1,5
1107                do j=1,nres_molec(i)
1108                itmp=itmp+1
1109               molnum(itmp)=i
1110                enddo
1111          enddo
1112 !        print *,nres_molec(i)
1113         endif
1114         print *,nres,"nres"
1115         if(.not.allocated(vbld)) then
1116            print *, "I DO ENTER"
1117            allocate(vbld(2*nres))
1118         endif
1119         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1120         do i=2,nres
1121           if (molnum(i).eq.1) then
1122           vbld(i)=vbl
1123           vbld_inv(i)=vblinv
1124
1125           else
1126           vbld(i)=7.0
1127           vbld_inv(i)=1.0/7.0
1128           endif
1129         enddo
1130         do i=2,nres-1
1131            if (molnum(i).eq.1) then
1132 !          print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i
1133           vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
1134           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
1135            else
1136           vbld(i+nres)=vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
1137           vbld_inv(i+nres)=1.0/vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
1138            endif
1139 !          write (iout,*) "i",i," itype",itype(i,1),
1140 !     &      " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1141         enddo
1142       endif
1143 !      print *,nres
1144 !      print '(20i4)',(itype(i,1),i=1,nres)
1145 !----------------------------
1146 !el reallocate tables
1147 !      do i=1,maxres2
1148 !        do j=1,3
1149 !          c_alloc(j,i)=c(j,i)
1150 !          dc_alloc(j,i)=dc(j,i)
1151 !        enddo
1152 !      enddo
1153 !      do i=1,maxres
1154 !elwrite(iout,*) "itype",i,itype(i,1)
1155 !        itype_alloc(i)=itype(i,1)
1156 !      enddo
1157
1158 !      deallocate(c)
1159 !      deallocate(dc)
1160 !      deallocate(itype)
1161 !      allocate(c(3,2*nres+4))
1162 !      allocate(dc(3,0:2*nres+2))
1163 !      allocate(itype(nres+2))
1164       allocate(itel(nres+2))
1165       itel(:)=0
1166
1167 !      do i=1,2*nres+2
1168 !        do j=1,3
1169 !          c(j,i)=c_alloc(j,i)
1170 !          dc(j,i)=dc_alloc(j,i)
1171 !        enddo
1172 !      enddo
1173 !      do i=1,nres+2
1174 !        itype(i,1)=itype_alloc(i)
1175 !        itel(i)=0
1176 !      enddo
1177 !--------------------------
1178       do i=1,nres
1179 #ifdef PROCOR
1180         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1181 #else
1182         if (itype(i,1).eq.ntyp1) then
1183 #endif
1184           itel(i)=0
1185 #ifdef PROCOR
1186         else if (iabs(itype(i+1,1)).ne.20) then
1187 #else
1188         else if (iabs(itype(i,1)).ne.20) then
1189 #endif
1190           itel(i)=1
1191         else
1192           itel(i)=2
1193         endif
1194       enddo
1195       if(me.eq.king.or..not.out1file)then
1196        write (iout,*) "ITEL"
1197        print *,nres,"nres"
1198        do i=1,nres-1
1199          write (iout,*) i,itype(i,1),itel(i)
1200        enddo
1201        print *,'Call Read_Bridge.'
1202       endif
1203       call read_bridge
1204 !--------------------------------
1205 !       print *,"tu dochodze"
1206 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1207       call alloc_geo_arrays
1208       call alloc_ener_arrays
1209 !--------------------------------
1210 ! 8/13/98 Set limits to generating the dihedral angles
1211       do i=1,nres
1212         phibound(1,i)=-pi
1213         phibound(2,i)=pi
1214       enddo
1215       read (inp,*) ndih_constr
1216       if (ndih_constr.gt.0) then
1217         raw_psipred=.false.
1218         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1219         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1220         allocate(ftors(ndih_constr)) !(maxdih_constr)
1221
1222 !        read (inp,*) ftors
1223         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), &
1224         i=1,ndih_constr)
1225         if(me.eq.king.or..not.out1file)then
1226          write (iout,*) &
1227          'There are',ndih_constr,' constraints on phi angles.'
1228          do i=1,ndih_constr
1229           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), &
1230           ftors(i)
1231          enddo
1232         endif
1233         do i=1,ndih_constr
1234           phi0(i)=deg2rad*phi0(i)
1235           drange(i)=deg2rad*drange(i)
1236         enddo
1237 !        if(me.eq.king.or..not.out1file) &
1238 !         write (iout,*) 'FTORS',ftors
1239         do i=1,ndih_constr
1240           ii = idih_constr(i)
1241           phibound(1,ii) = phi0(i)-drange(i)
1242           phibound(2,ii) = phi0(i)+drange(i)
1243         enddo
1244       else if (ndih_constr.lt.0) then
1245         raw_psipred=.true.
1246         allocate(idih_constr(nres))
1247         allocate(secprob(3,nres))
1248         allocate(vpsipred(3,nres))
1249         allocate(sdihed(2,nres))
1250         call card_concat(weightcard,.true.)
1251         call reada(weightcard,"PHIHEL",phihel,50.0D0)
1252         call reada(weightcard,"PHIBET",phibet,180.0D0)
1253         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
1254         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
1255         call reada(weightcard,"WDIHC",wdihc,0.591D0)
1256         write (iout,*) "Weight of dihedral angle restraints",wdihc
1257         read(inp,'(9x,3f7.3)') &
1258           (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
1259         write (iout,*) "The secprob array"
1260         do i=nnt,nct
1261           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
1262         enddo
1263         ndih_constr=0
1264         do i=nnt+3,nct
1265           if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 &
1266           .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then
1267             ndih_constr=ndih_constr+1
1268             idih_constr(ndih_constr)=i
1269             sumv=0.0d0
1270             do j=1,3
1271               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
1272               sumv=sumv+vpsipred(j,ndih_constr)
1273             enddo
1274             do j=1,3
1275               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1276             enddo
1277             phibound(1,ndih_constr)=phihel*deg2rad
1278             phibound(2,ndih_constr)=phibet*deg2rad
1279             sdihed(1,ndih_constr)=sigmahel*deg2rad
1280             sdihed(2,ndih_constr)=sigmabet*deg2rad
1281           endif
1282         enddo
1283
1284       endif
1285       if (with_theta_constr) then
1286 !C with_theta_constr is keyword allowing for occurance of theta constrains
1287       read (inp,*) ntheta_constr
1288 !C ntheta_constr is the number of theta constrains
1289       if (ntheta_constr.gt.0) then
1290 !C        read (inp,*) ftors
1291         allocate(itheta_constr(ntheta_constr))
1292         allocate(theta_constr0(ntheta_constr))
1293         allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1294         read (inp,*) (itheta_constr(i),theta_constr0(i), &
1295        theta_drange(i),for_thet_constr(i), &
1296        i=1,ntheta_constr)
1297 !C the above code reads from 1 to ntheta_constr
1298 !C itheta_constr(i) residue i for which is theta_constr
1299 !C theta_constr0 the global minimum value
1300 !C theta_drange is range for which there is no energy penalty
1301 !C for_thet_constr is the force constant for quartic energy penalty
1302 !C E=k*x**4
1303         if(me.eq.king.or..not.out1file)then
1304          write (iout,*) &
1305         'There are',ntheta_constr,' constraints on phi angles.'
1306          do i=1,ntheta_constr
1307           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1308          theta_drange(i), &
1309          for_thet_constr(i)
1310          enddo
1311         endif
1312         do i=1,ntheta_constr
1313           theta_constr0(i)=deg2rad*theta_constr0(i)
1314           theta_drange(i)=deg2rad*theta_drange(i)
1315         enddo
1316 !C        if(me.eq.king.or..not.out1file)
1317 !C     &   write (iout,*) 'FTORS',ftors
1318 !C        do i=1,ntheta_constr
1319 !C          ii = itheta_constr(i)
1320 !C          thetabound(1,ii) = phi0(i)-drange(i)
1321 !C          thetabound(2,ii) = phi0(i)+drange(i)
1322 !C        enddo
1323       endif ! ntheta_constr.gt.0
1324       endif! with_theta_constr
1325
1326       nnt=1
1327 #ifdef MPI
1328       if (me.eq.king) then
1329 #endif
1330        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1331        do i=1,nres
1332          write (iout,'(a3,i5,2f10.1)') &
1333          restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1334        enddo
1335 #ifdef MP
1336       endif
1337 #endif
1338       nct=nres
1339       print *,'NNT=',NNT,' NCT=',NCT
1340       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1341       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1342       if (pdbref) then
1343         if(me.eq.king.or..not.out1file) &
1344          write (iout,'(a,i3)') 'nsup=',nsup
1345         nstart_seq=nnt
1346         if (nsup.le.(nct-nnt+1)) then
1347           do i=0,nct-nnt+1-nsup
1348             if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1349               nstart_seq=nnt+i
1350               goto 111
1351             endif
1352           enddo
1353           write (iout,'(a)') &
1354                   'Error - sequences to be superposed do not match.'
1355           stop
1356         else
1357           do i=0,nsup-(nct-nnt+1)
1358             if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1359             then
1360               nstart_sup=nstart_sup+i
1361               nsup=nct-nnt+1
1362               goto 111
1363             endif
1364           enddo
1365           write (iout,'(a)') &
1366                   'Error - sequences to be superposed do not match.'
1367         endif
1368   111   continue
1369         if (nsup.eq.0) nsup=nct-nnt
1370         if (nstart_sup.eq.0) nstart_sup=nnt
1371         if (nstart_seq.eq.0) nstart_seq=nnt
1372         if(me.eq.king.or..not.out1file) &
1373          write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1374                        ' nstart_seq=',nstart_seq !,"242343453254"
1375       endif
1376 !--- Zscore rms -------
1377       if (nz_start.eq.0) nz_start=nnt
1378       if (nz_end.eq.0 .and. nsup.gt.0) then
1379         nz_end=nnt+nsup-1
1380       else if (nz_end.eq.0) then
1381         nz_end=nct
1382       endif
1383       if(me.eq.king.or..not.out1file)then
1384        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1385        write (iout,*) 'IZ_SC=',iz_sc
1386       endif
1387 !----------------------
1388       call init_int_table
1389       if (refstr) then
1390         if (.not.pdbref) then
1391           call read_angles(inp,*38)
1392           goto 39
1393    38     write (iout,'(a)') 'Error reading reference structure.'
1394 #ifdef MPI
1395           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1396           stop 'Error reading reference structure'
1397 #endif
1398    39     call chainbuild
1399           call setup_var
1400 !zscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1401           nstart_sup=nnt
1402           nstart_seq=nnt
1403           nsup=nct-nnt+1
1404           kkk=1
1405           do i=1,2*nres
1406             do j=1,3
1407               cref(j,i,kkk)=c(j,i)
1408             enddo
1409           enddo
1410           call contact(.true.,ncont_ref,icont_ref,co)
1411         endif
1412 !        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1413 !        call flush(iout)
1414 !EL        if (constr_dist.gt.0) call read_dist_constr
1415 !EL        write (iout,*) "After read_dist_constr nhpb",nhpb
1416 !EL        if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1417 !EL        call hpb_partition
1418         if(me.eq.king.or..not.out1file) &
1419          write (iout,*) 'Contact order:',co
1420         if (pdbref) then
1421         if(me.eq.king.or..not.out1file) &
1422          write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1423         do i=1,ncont_ref
1424           do j=1,2
1425             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1426           enddo
1427           if(me.eq.king.or..not.out1file) &
1428            write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1429            icont_ref(1,i),' ',&
1430            restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1431         enddo
1432         endif
1433       if (constr_homology.gt.0) then
1434 !        write (iout,*) "Calling read_constr_homology"
1435 !        call flush(iout)
1436         call read_constr_homology
1437         if (indpdb.gt.0 .or. pdbref) then
1438           do i=1,2*nres
1439             do j=1,3
1440               c(j,i)=crefjlee(j,i)
1441               cref(j,i,1)=crefjlee(j,i)
1442             enddo
1443           enddo
1444         endif
1445 #define DEBUG
1446 #ifdef DEBUG
1447         write (iout,*) "sc_loc_geom: Array C"
1448         do i=1,nres
1449           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
1450            (c(j,i+nres),j=1,3)
1451         enddo
1452         write (iout,*) "Array Cref"
1453         do i=1,nres
1454           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),&
1455            (cref(j,i+nres,1),j=1,3)
1456         enddo
1457 #endif
1458        call int_from_cart1(.false.)
1459        call sc_loc_geom(.false.)
1460        do i=1,nres
1461          thetaref(i)=theta(i)
1462          phiref(i)=phi(i)
1463        enddo
1464        do i=1,nres-1
1465          do j=1,3
1466            dc(j,i)=c(j,i+1)-c(j,i)
1467            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1468          enddo
1469        enddo
1470        do i=2,nres-1
1471          do j=1,3
1472            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1473            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1474          enddo
1475        enddo
1476       else
1477         homol_nset=0
1478         if (start_from_model) then
1479           nmodel_start=0
1480           do
1481             read(inp,'(a)',end=332,err=332) pdbfile
1482             if (me.eq.king .or. .not. out1file)&
1483              write (iout,'(a,5x,a)') 'Opening PDB file',&
1484              pdbfile(:ilen(pdbfile))
1485             open(ipdbin,file=pdbfile,status='old',err=336)
1486             goto 335
1487  336        write (iout,'(a,5x,a)') 'Error opening PDB file',&
1488            pdbfile(:ilen(pdbfile))
1489             call flush(iout)
1490             stop
1491  335        continue
1492             unres_pdb=.false.
1493             nres_temp=nres
1494 !            call readpdb
1495             call readpdb_template(nmodel_start+1)
1496             close(ipdbin)
1497             if (nres.ge.nres_temp) then
1498               nmodel_start=nmodel_start+1
1499               pdbfiles_chomo(nmodel_start)=pdbfile
1500               do i=1,2*nres
1501                 do j=1,3
1502                   chomo(j,i,nmodel_start)=c(j,i)
1503                 enddo
1504               enddo
1505             else
1506               if (me.eq.king .or. .not. out1file) &
1507                write (iout,'(a,2i7,1x,a)') &
1508                 "Different number of residues",nres_temp,nres, &
1509                 " model skipped."
1510             endif
1511             nres=nres_temp
1512           enddo
1513   332     continue
1514           if (nmodel_start.eq.0) then
1515             if (me.eq.king .or. .not. out1file) &
1516              write (iout,'(a)') &
1517              "No valid starting model found START_FROM_MODELS is OFF"
1518               start_from_model=.false.
1519           endif
1520           write (iout,*) "nmodel_start",nmodel_start
1521         endif
1522       endif
1523
1524       endif
1525         if (constr_dist.gt.0) call read_dist_constr
1526         write (iout,*) "After read_dist_constr nhpb",nhpb
1527         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1528         call hpb_partition
1529
1530       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1531           .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1532           modecalc.ne.10) then
1533 ! If input structure hasn't been supplied from the PDB file read or generate
1534 ! initial geometry.
1535         if (iranconf.eq.0 .and. .not. extconf) then
1536           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1537            write (iout,'(a)') 'Initial geometry will be read in.'
1538           if (read_cart) then
1539             read(inp,'(8f10.5)',end=36,err=36) &
1540              ((c(l,k),l=1,3),k=1,nres),&
1541              ((c(l,k+nres),l=1,3),k=nnt,nct)
1542             write (iout,*) "Exit READ_CART"
1543             write (iout,'(8f10.5)') &
1544              ((c(l,k),l=1,3),k=1,nres)
1545             write (iout,'(8f10.5)') &
1546              ((c(l,k+nres),l=1,3),k=nnt,nct)
1547             call int_from_cart1(.true.)
1548             write (iout,*) "Finish INT_TO_CART"
1549             do i=1,nres-1
1550               do j=1,3
1551                 dc(j,i)=c(j,i+1)-c(j,i)
1552                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1553               enddo
1554             enddo
1555             do i=nnt,nct
1556               if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1557                 do j=1,3
1558                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1559                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1560                 enddo
1561               endif
1562             enddo
1563             return
1564           else
1565            write(iout,*) "read angles from input"
1566            call read_angles(inp,*36)
1567             call chainbuild
1568
1569           endif
1570           goto 37
1571    36     write (iout,'(a)') 'Error reading angle file.'
1572 #ifdef MPI
1573           call mpi_finalize( MPI_COMM_WORLD,IERR )
1574 #endif
1575           stop 'Error reading angle file.'
1576    37     continue
1577         else if (extconf) then
1578          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1579           write (iout,'(a)') 'Extended chain initial geometry.'
1580          do i=3,nres
1581           theta(i)=90d0*deg2rad
1582          enddo
1583          do i=4,nres
1584           phi(i)=180d0*deg2rad
1585          enddo
1586          do i=2,nres-1
1587           alph(i)=110d0*deg2rad
1588          enddo
1589          do i=2,nres-1
1590           omeg(i)=-120d0*deg2rad
1591           if (itype(i,1).le.0) omeg(i)=-omeg(i)
1592          enddo
1593          call chainbuild
1594         else
1595           if(me.eq.king.or..not.out1file) &
1596            write (iout,'(a)') 'Random-generated initial geometry.'
1597
1598
1599 #ifdef MPI
1600           if (me.eq.king  .or. fg_rank.eq.0 .and. &
1601                  ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1602 #endif
1603             do itrial=1,100
1604               itmp=1
1605               call gen_rand_conf(itmp,*30)
1606               goto 40
1607    30         write (iout,*) 'Failed to generate random conformation',&
1608                 ', itrial=',itrial
1609               write (*,*) 'Processor:',me,&
1610                 ' Failed to generate random conformation',&
1611                 ' itrial=',itrial
1612               call intout
1613
1614 #ifdef AIX
1615               call flush_(iout)
1616 #else
1617               call flush(iout)
1618 #endif
1619             enddo
1620             write (iout,'(a,i3,a)') 'Processor:',me,&
1621               ' error in generating random conformation.'
1622             write (*,'(a,i3,a)') 'Processor:',me,&
1623               ' error in generating random conformation.'
1624             call flush(iout)
1625 #ifdef MPI
1626             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1627    40       continue
1628           endif
1629 #else
1630           do itrial=1,100
1631             itmp=1
1632             call gen_rand_conf(itmp,*335)
1633             goto 40
1634   335       write (iout,*) 'Failed to generate random conformation',&
1635               ', itrial=',itrial
1636             write (*,*) 'Failed to generate random conformation',&
1637               ', itrial=',itrial
1638           enddo
1639           write (iout,'(a,i3,a)') 'Processor:',me,&
1640             ' error in generating random conformation.'
1641           write (*,'(a,i3,a)') 'Processor:',me,&
1642             ' error in generating random conformation.'
1643           stop
1644    40     continue
1645 #endif
1646         endif
1647       elseif (modecalc.eq.4) then
1648         read (inp,'(a)') intinname
1649         open (intin,file=intinname,status='old',err=333)
1650         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1651         write (iout,'(a)') 'intinname',intinname
1652         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1653         goto 334
1654   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1655 #ifdef MPI
1656         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1657 #endif
1658         stop 'Error opening angle file.'
1659   334   continue
1660
1661       endif
1662 ! Generate distance constraints, if the PDB structure is to be regularized.
1663       if (nthread.gt.0) then
1664         call read_threadbase
1665       endif
1666       call setup_var
1667       if (me.eq.king .or. .not. out1file) &
1668        call intout
1669       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1670         write (iout,'(/a,i3,a)') &
1671         'The chain contains',ns,' disulfide-bridging cysteines.'
1672         write (iout,'(20i4)') (iss(i),i=1,ns)
1673        if (dyn_ss) then
1674           write(iout,*)"Running with dynamic disulfide-bond formation"
1675        else
1676         write (iout,'(/a/)') 'Pre-formed links are:'
1677         do i=1,nss
1678           i1=ihpb(i)-nres
1679           i2=jhpb(i)-nres
1680           it1=itype(i1,1)
1681           it2=itype(i2,1)
1682           if (me.eq.king.or..not.out1file) &
1683           write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1684           restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1685           ebr,forcon(i)
1686         enddo
1687         write (iout,'(a)')
1688        endif
1689       endif
1690       if (ns.gt.0.and.dyn_ss) then
1691           do i=nss+1,nhpb
1692             ihpb(i-nss)=ihpb(i)
1693             jhpb(i-nss)=jhpb(i)
1694             forcon(i-nss)=forcon(i)
1695             dhpb(i-nss)=dhpb(i)
1696           enddo
1697           nhpb=nhpb-nss
1698           nss=0
1699           call hpb_partition
1700           do i=1,ns
1701             dyn_ss_mask(iss(i))=.true.
1702           enddo
1703       endif
1704       if (i2ndstr.gt.0) call secstrp2dihc
1705       if (indpdb.gt.0) then
1706           write(iout,*) "WCHODZE TU!!"
1707           call int_from_cart1(.true.)
1708       endif
1709 !      call geom_to_var(nvar,x)
1710 !      call etotal(energia(0))
1711 !      call enerprint(energia(0))
1712 !      call briefout(0,etot)
1713 !      stop
1714 !d    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1715 !d    write (iout,'(a)') 'Variable list:'
1716 !d    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1717 #ifdef MPI
1718       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1719         write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1720         'Processor',myrank,': end reading molecular data.'
1721 #endif
1722       return
1723       end subroutine molread
1724 !-----------------------------------------------------------------------------
1725       subroutine read_constr_homology
1726       use control, only:init_int_table,homology_partition
1727       use MD_data, only:iset
1728 !      implicit none
1729 !      include 'DIMENSIONS'
1730 !#ifdef MPI
1731 !      include 'mpif.h'
1732 !#endif
1733 !      include 'COMMON.SETUP'
1734 !      include 'COMMON.CONTROL'
1735 !      include 'COMMON.HOMOLOGY'
1736 !      include 'COMMON.CHAIN'
1737 !      include 'COMMON.IOUNITS'
1738 !      include 'COMMON.MD'
1739 !      include 'COMMON.QRESTR'
1740 !      include 'COMMON.GEO'
1741 !      include 'COMMON.INTERACT'
1742 !      include 'COMMON.NAMES'
1743 !      include 'COMMON.VAR'
1744 !
1745
1746 !     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1747 !    &                 dist_cut
1748 !     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1749 !    &    sigma_odl_temp(maxres,maxres,max_template)
1750       character*2 kic2
1751       character*24 model_ki_dist, model_ki_angle
1752       character*500 controlcard
1753       integer :: ki,i,ii,j,k,l
1754       integer, dimension (:), allocatable :: ii_in_use
1755       integer :: i_tmp,idomain_tmp,&
1756       irec,ik,iistart,nres_temp
1757 !      integer :: iset
1758 !      external :: ilen
1759       logical :: liiflag,lfirst
1760       integer :: i01,i10
1761 !
1762 !     FP - Nov. 2014 Temporary specifications for new vars
1763 !
1764       real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1765                        rescore3_tmp, dist_cut
1766       real(kind=8), dimension (:,:),allocatable :: rescore
1767       real(kind=8), dimension (:,:),allocatable :: rescore2
1768       real(kind=8), dimension (:,:),allocatable :: rescore3
1769       real(kind=8) :: distal
1770       character*24 tpl_k_rescore
1771       character*256 pdbfile
1772
1773 ! -----------------------------------------------------------------
1774 ! Reading multiple PDB ref structures and calculation of retraints
1775 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1776 ! FP (Nov., 2014)
1777 ! -----------------------------------------------------------------
1778 !
1779 !
1780 ! Alternative: reading from input
1781       call card_concat(controlcard,.true.)
1782       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1783       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1784       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1785       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1786       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1787       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1788       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1789       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1790       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1791       if(.not.read2sigma.and.start_from_model) then
1792           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1793            write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1794           start_from_model=.false.
1795           iranconf=(indpdb.le.0)
1796       endif
1797       if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1798          write(iout,*) 'START_FROM_MODELS is ON'
1799 !      if(start_from_model .and. rest) then
1800 !        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1801 !         write(iout,*) 'START_FROM_MODELS is OFF'
1802 !         write(iout,*) 'remove restart keyword from input'
1803 !        endif
1804 !      endif
1805       if (start_from_model) nmodel_start=constr_homology
1806       if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1807       if (homol_nset.gt.1)then
1808          call card_concat(controlcard,.true.)
1809          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1810          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1811 !          write(iout,*) "iset homology_weight "
1812           do i=1,homol_nset
1813            write(iout,*) i,waga_homology(i)
1814           enddo
1815          endif
1816          iset=mod(kolor,homol_nset)+1
1817       else
1818        iset=1
1819        waga_homology(1)=1.0
1820       endif
1821
1822 !d      write (iout,*) "nnt",nnt," nct",nct
1823 !d      call flush(iout)
1824
1825       if (read_homol_frag) then
1826         call read_klapaucjusz
1827       else
1828
1829       lim_odl=0
1830       lim_dih=0
1831 !
1832 !      write(iout,*) 'nnt=',nnt,'nct=',nct
1833 !
1834 !      do i = nnt,nct
1835 !        do k=1,constr_homology
1836 !         idomain(k,i)=0
1837 !        enddo
1838 !      enddo
1839        idomain=0
1840
1841 !      ii=0
1842 !      do i = nnt,nct-2
1843 !        do j=i+2,nct
1844 !        ii=ii+1
1845 !        ii_in_use(ii)=0
1846 !        enddo
1847 !      enddo
1848       ii_in_use=0
1849       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1850       if(.not.allocated(chomo)) allocate(chomo(3,nres,constr_homology))
1851
1852       do k=1,constr_homology
1853
1854         read(inp,'(a)') pdbfile
1855         pdbfiles_chomo(k)=pdbfile
1856         if(me.eq.king .or. .not. out1file) &
1857          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1858         pdbfile(:ilen(pdbfile))
1859         open(ipdbin,file=pdbfile,status='old',err=33)
1860         goto 34
1861   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
1862         pdbfile(:ilen(pdbfile))
1863         stop
1864   34    continue
1865 !        print *,'Begin reading pdb data'
1866 !
1867 ! Files containing res sim or local scores (former containing sigmas)
1868 !
1869
1870         write(kic2,'(bz,i2.2)') k
1871
1872         tpl_k_rescore="template"//kic2//".sco"
1873         write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1874         unres_pdb=.false.
1875         nres_temp=nres
1876         write(iout,*) "read2sigma",read2sigma
1877
1878         if (read2sigma) then
1879           call readpdb_template(k)
1880         else
1881           call readpdb
1882         endif
1883         write(iout,*) "after readpdb"
1884         if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1885         nres_chomo(k)=nres
1886         nres=nres_temp
1887         if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1888         if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1889         if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1890         if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1891         if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1892         if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1893         if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1894         if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1895         if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1896         if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1897         if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1898         if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1899         if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1900         if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1901 !        if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1902         if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1903         if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1904         if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1905         if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1906 !        if(.not.allocated(distance)) allocate(distance(constr_homology))
1907 !        if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1908
1909
1910 !
1911 !     Distance restraints
1912 !
1913 !          ... --> odl(k,ii)
1914 ! Copy the coordinates from reference coordinates (?)
1915         do i=1,2*nres_chomo(k)
1916           do j=1,3
1917             c(j,i)=cref(j,i,1)
1918 !           write (iout,*) "c(",j,i,") =",c(j,i)
1919           enddo
1920         enddo
1921 !
1922 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1923 !
1924 !         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1925           open (ientin,file=tpl_k_rescore,status='old')
1926           if (nnt.gt.1) rescore(k,1)=0.0d0
1927           do irec=nnt,nct ! loop for reading res sim
1928             if (read2sigma) then
1929              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1930                                      rescore3_tmp,idomain_tmp
1931              i_tmp=i_tmp+nnt-1
1932              idomain(k,i_tmp)=idomain_tmp
1933              rescore(k,i_tmp)=rescore_tmp
1934              rescore2(k,i_tmp)=rescore2_tmp
1935              rescore3(k,i_tmp)=rescore3_tmp
1936              if (.not. out1file .or. me.eq.king)&
1937              write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
1938                            i_tmp,rescore2_tmp,rescore_tmp,&
1939                                      rescore3_tmp,idomain_tmp
1940             else
1941              idomain(k,irec)=1
1942              read (ientin,*,end=1401) rescore_tmp
1943
1944 !           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1945              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1946 !           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1947             endif
1948           enddo
1949  1401   continue
1950         close (ientin)
1951         if (waga_dist.ne.0.0d0) then
1952           ii=0
1953           do i = nnt,nct-2
1954             do j=i+2,nct
1955
1956               x12=c(1,i)-c(1,j)
1957               y12=c(2,i)-c(2,j)
1958               z12=c(3,i)-c(3,j)
1959               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1960 !              write (iout,*) k,i,j,distal,dist2_cut
1961
1962             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
1963                  .and. distal.le.dist2_cut ) then
1964
1965               ii=ii+1
1966               ii_in_use(ii)=1
1967               l_homo(k,ii)=.true.
1968
1969 !             write (iout,*) "k",k
1970 !             write (iout,*) "i",i," j",j," constr_homology",
1971 !    &                       constr_homology
1972               ires_homo(ii)=i
1973               jres_homo(ii)=j
1974               odl(k,ii)=distal
1975               if (read2sigma) then
1976                 sigma_odl(k,ii)=0
1977                 do ik=i,j
1978                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1979                 enddo
1980                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1981                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
1982               sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1983               else
1984                 if (odl(k,ii).le.dist_cut) then
1985                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1986                 else
1987 #ifdef OLDSIGMA
1988                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
1989                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1990 #else
1991                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
1992                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1993 #endif
1994                 endif
1995               endif
1996               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1997             else
1998 !              ii=ii+1
1999 !              l_homo(k,ii)=.false.
2000             endif
2001             enddo
2002           enddo
2003         lim_odl=ii
2004         endif
2005 !        write (iout,*) "Distance restraints set"
2006 !        call flush(iout)
2007 !
2008 !     Theta, dihedral and SC retraints
2009 !
2010         if (waga_angle.gt.0.0d0) then
2011 !         open (ientin,file=tpl_k_sigma_dih,status='old')
2012 !         do irec=1,maxres-3 ! loop for reading sigma_dih
2013 !            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2014 !            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2015 !            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2016 !    &                            sigma_dih(k,i+nnt-1)
2017 !         enddo
2018 !1402   continue
2019 !         close (ientin)
2020           do i = nnt+3,nct
2021             if (idomain(k,i).eq.0) then
2022                sigma_dih(k,i)=0.0
2023                cycle
2024             endif
2025             dih(k,i)=phiref(i) ! right?
2026 !           read (ientin,*) sigma_dih(k,i) ! original variant
2027 !             write (iout,*) "dih(",k,i,") =",dih(k,i)
2028 !             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2029 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2030 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
2031 !    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
2032
2033             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2034                           rescore(k,i-2)+rescore(k,i-3))/4.0
2035 !            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2036 !           write (iout,*) "Raw sigmas for dihedral angle restraints"
2037 !           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2038 !           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2039 !                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
2040 !   Instead of res sim other local measure of b/b str reliability possible
2041             if (sigma_dih(k,i).ne.0) &
2042             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2043 !           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2044           enddo
2045           lim_dih=nct-nnt-2
2046         endif
2047 !        write (iout,*) "Dihedral angle restraints set"
2048 !        call flush(iout)
2049
2050         if (waga_theta.gt.0.0d0) then
2051 !         open (ientin,file=tpl_k_sigma_theta,status='old')
2052 !         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2053 !            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2054 !            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2055 !    &                              sigma_theta(k,i+nnt-1)
2056 !         enddo
2057 !1403   continue
2058 !         close (ientin)
2059
2060           do i = nnt+2,nct ! right? without parallel.
2061 !         do i = i=1,nres ! alternative for bounds acc to readpdb?
2062 !         do i=ithet_start,ithet_end ! with FG parallel.
2063              if (idomain(k,i).eq.0) then
2064               sigma_theta(k,i)=0.0
2065               cycle
2066              endif
2067              thetatpl(k,i)=thetaref(i)
2068 !            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2069 !            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
2070 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2071 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
2072 !            read (ientin,*) sigma_theta(k,i) ! 1st variant
2073              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2074                              rescore(k,i-2))/3.0
2075 !             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2076              if (sigma_theta(k,i).ne.0) &
2077              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2078
2079 !            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2080 !                             rescore(k,i-2) !  right expression ?
2081 !            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2082           enddo
2083         endif
2084 !        write (iout,*) "Angle restraints set"
2085 !        call flush(iout)
2086
2087         if (waga_d.gt.0.0d0) then
2088 !       open (ientin,file=tpl_k_sigma_d,status='old')
2089 !         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2090 !            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2091 !            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2092 !    &                          sigma_d(k,i+nnt-1)
2093 !         enddo
2094 !1404   continue
2095
2096           do i = nnt,nct ! right? without parallel.
2097 !         do i=2,nres-1 ! alternative for bounds acc to readpdb?
2098 !         do i=loc_start,loc_end ! with FG parallel.
2099                if (itype(i,1).eq.10) cycle
2100                if (idomain(k,i).eq.0 ) then
2101                   sigma_d(k,i)=0.0
2102                   cycle
2103                endif
2104                xxtpl(k,i)=xxref(i)
2105                yytpl(k,i)=yyref(i)
2106                zztpl(k,i)=zzref(i)
2107 !              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2108 !              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2109 !              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2110 !              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
2111                sigma_d(k,i)=rescore3(k,i) !  right expression ?
2112                if (sigma_d(k,i).ne.0) &
2113                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2114
2115 !              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
2116 !              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2117 !              read (ientin,*) sigma_d(k,i) ! 1st variant
2118           enddo
2119         endif
2120       enddo
2121 !      write (iout,*) "SC restraints set"
2122 !      call flush(iout)
2123 !
2124 ! remove distance restraints not used in any model from the list
2125 ! shift data in all arrays
2126 !
2127 !      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
2128       if (waga_dist.ne.0.0d0) then
2129         ii=0
2130         liiflag=.true.
2131         lfirst=.true.
2132         do i=nnt,nct-2
2133          do j=i+2,nct
2134           ii=ii+1
2135 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2136 !     &            .and. distal.le.dist2_cut ) then
2137 !          write (iout,*) "i",i," j",j," ii",ii
2138 !          call flush(iout)
2139           if (ii_in_use(ii).eq.0.and.liiflag.or. &
2140           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2141             liiflag=.false.
2142             i10=ii
2143             if (lfirst) then
2144               lfirst=.false.
2145               iistart=ii
2146             else
2147               if(i10.eq.lim_odl) i10=i10+1
2148               do ki=0,i10-i01-1
2149                ires_homo(iistart+ki)=ires_homo(ki+i01)
2150                jres_homo(iistart+ki)=jres_homo(ki+i01)
2151                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2152                do k=1,constr_homology
2153                 odl(k,iistart+ki)=odl(k,ki+i01)
2154                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2155                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2156                enddo
2157               enddo
2158               iistart=iistart+i10-i01
2159             endif
2160           endif
2161           if (ii_in_use(ii).ne.0.and..not.liiflag) then
2162              i01=ii
2163              liiflag=.true.
2164           endif
2165          enddo
2166         enddo
2167         lim_odl=iistart-1
2168       endif
2169 !      write (iout,*) "Removing distances completed"
2170 !      call flush(iout)
2171       endif ! .not. klapaucjusz
2172
2173       if (constr_homology.gt.0) call homology_partition
2174       write (iout,*) "After homology_partition"
2175       call flush(iout)
2176       if (constr_homology.gt.0) call init_int_table
2177       write (iout,*) "After init_int_table"
2178       call flush(iout)
2179 !      endif ! .not. klapaucjusz
2180 !      endif
2181 !      if (constr_homology.gt.0) call homology_partition
2182 !      write (iout,*) "After homology_partition"
2183 !      call flush(iout)
2184 !      if (constr_homology.gt.0) call init_int_table
2185 !      write (iout,*) "After init_int_table"
2186 !      call flush(iout)
2187 !      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2188 !      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2189 !
2190 ! Print restraints
2191 !
2192       if (.not.out_template_restr) return
2193 !d      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2194       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2195        write (iout,*) "Distance restraints from templates"
2196        do ii=1,lim_odl
2197        write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
2198         ii,ires_homo(ii),jres_homo(ii),&
2199         (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
2200         ki=1,constr_homology)
2201        enddo
2202        write (iout,*) "Dihedral angle restraints from templates"
2203        do i=nnt+3,nct
2204         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2205             (rad2deg*dih(ki,i),&
2206             rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2207        enddo
2208        write (iout,*) "Virtual-bond angle restraints from templates"
2209        do i=nnt+2,nct
2210         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2211             (rad2deg*thetatpl(ki,i),&
2212             rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2213        enddo
2214        write (iout,*) "SC restraints from templates"
2215        do i=nnt,nct
2216         write(iout,'(i7,100(4f8.2,4x))') i,&
2217         (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
2218          1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2219        enddo
2220       endif
2221       return
2222       end subroutine read_constr_homology
2223 !-----------------------------------------------------------------------------
2224       subroutine read_klapaucjusz
2225 !     implicit none
2226 !     include 'DIMENSIONS'
2227 !#ifdef MPI
2228 !     include 'mpif.h'
2229 !#endif
2230 !     include 'COMMON.SETUP'
2231 !     include 'COMMON.CONTROL'
2232 !     include 'COMMON.HOMOLOGY'
2233 !     include 'COMMON.CHAIN'
2234 !     include 'COMMON.IOUNITS'
2235 !     include 'COMMON.MD'
2236 !     include 'COMMON.GEO'
2237 !     include 'COMMON.INTERACT'
2238 !     include 'COMMON.NAMES'
2239       character(len=256) fragfile
2240       integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2241                          ii_in_use
2242       integer, dimension(:,:), allocatable :: iresclust,inclust
2243       integer :: nclust
2244
2245       character(len=2) :: kic2
2246       character(len=24) :: model_ki_dist, model_ki_angle
2247       character(len=500) :: controlcard
2248       integer :: ki, i, j, jj,k, l, i_tmp,&
2249       idomain_tmp,&
2250       ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2251       i01,i10,nnt_chain,nct_chain
2252       real(kind=8) :: distal
2253       logical :: lprn = .true.
2254       integer :: nres_temp
2255 !      integer :: ilen
2256 !      external :: ilen
2257       logical :: liiflag,lfirst
2258
2259       real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2260       real(kind=8), dimension (:,:), allocatable  :: rescore
2261       real(kind=8), dimension (:,:), allocatable :: rescore2
2262       character(len=24) :: tpl_k_rescore
2263       character(len=256) :: pdbfile
2264
2265 !
2266 ! For new homol impl
2267 !
2268 !     include 'COMMON.VAR'
2269 !
2270 !      write (iout,*) "READ_KLAPAUCJUSZ"
2271 !      print *,"READ_KLAPAUCJUSZ"
2272 !      call flush(iout)
2273       call getenv("FRAGFILE",fragfile)
2274       write (iout,*) "Opening", fragfile
2275       call flush(iout)
2276       open(ientin,file=fragfile,status="old",err=10)
2277 !      write (iout,*) " opened"
2278 !      call flush(iout)
2279
2280       sigma_theta=0.0
2281       sigma_d=0.0
2282       sigma_dih=0.0
2283       l_homo = .false.
2284
2285       nres_temp=nres
2286       itype_temp(:)=itype(:,1)
2287       ii=0
2288       lim_odl=0
2289
2290 !      write (iout,*) "Entering loop"
2291 !      call flush(iout)
2292
2293       DO IGR = 1,NCHAIN_GROUP
2294
2295 !      write (iout,*) "igr",igr
2296       call flush(iout)
2297       read(ientin,*) constr_homology,nclust
2298       if (start_from_model) then
2299         nmodel_start=constr_homology
2300       else
2301         nmodel_start=0
2302       endif
2303
2304       ii_old=lim_odl
2305
2306       ichain=iequiv(1,igr)
2307       nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2308       nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2309 !      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2310
2311 ! Read pdb files
2312       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2313       if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2314       do k=1,constr_homology
2315         read(ientin,'(a)') pdbfile
2316         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2317         pdbfile(:ilen(pdbfile))
2318         pdbfiles_chomo(k)=pdbfile
2319         open(ipdbin,file=pdbfile,status='old',err=33)
2320         goto 34
2321   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
2322         pdbfile(:ilen(pdbfile))
2323         stop
2324   34    continue
2325         unres_pdb=.false.
2326 !        nres_temp=nres
2327         call readpdb_template(k)
2328         nres_chomo(k)=nres
2329 !        nres=nres_temp
2330         do i=1,nres
2331           rescore(k,i)=0.2d0
2332           rescore2(k,i)=1.0d0
2333         enddo
2334       enddo
2335 ! Read clusters
2336       do i=1,nclust
2337         read(ientin,*) ninclust(i),nresclust(i)
2338         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2339         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2340       enddo
2341 !
2342 ! Loop over clusters
2343 !
2344       do l=1,nclust
2345         do ll = 1,ninclust(l)
2346
2347         k = inclust(ll,l)
2348 !        write (iout,*) "l",l," ll",ll," k",k
2349         do i=1,nres
2350           idomain(k,i)=0
2351         enddo
2352         do i=1,nresclust(l)
2353           if (nnt.gt.1)  then
2354             idomain(k,iresclust(i,l)+1) = 1
2355           else
2356             idomain(k,iresclust(i,l)) = 1
2357           endif
2358         enddo
2359 !
2360 !     Distance restraints
2361 !
2362 !          ... --> odl(k,ii)
2363 ! Copy the coordinates from reference coordinates (?)
2364 !        nres_temp=nres
2365         nres=nres_chomo(k)
2366         do i=1,2*nres
2367           do j=1,3
2368             c(j,i)=chomo(j,i,k)
2369 !           write (iout,*) "c(",j,i,") =",c(j,i)
2370           enddo
2371         enddo
2372         call int_from_cart(.true.,.false.)
2373         call sc_loc_geom(.false.)
2374         do i=1,nres
2375           thetaref(i)=theta(i)
2376           phiref(i)=phi(i)
2377         enddo
2378 !        nres=nres_temp
2379         if (waga_dist.ne.0.0d0) then
2380           ii=ii_old
2381 !          do i = nnt,nct-2
2382           do i = nnt_chain,nct_chain-2
2383 !            do j=i+2,nct
2384             do j=i+2,nct_chain
2385
2386               x12=c(1,i)-c(1,j)
2387               y12=c(2,i)-c(2,j)
2388               z12=c(3,i)-c(3,j)
2389               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2390 !              write (iout,*) k,i,j,distal,dist2_cut
2391
2392             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2393                  .and. distal.le.dist2_cut ) then
2394
2395               ii=ii+1
2396               ii_in_use(ii)=1
2397               l_homo(k,ii)=.true.
2398
2399 !             write (iout,*) "k",k
2400 !             write (iout,*) "i",i," j",j," constr_homology",
2401 !     &                       constr_homology
2402               ires_homo(ii)=i+chain_border1(1,igr)-1
2403               jres_homo(ii)=j+chain_border1(1,igr)-1
2404               odl(k,ii)=distal
2405               if (read2sigma) then
2406                 sigma_odl(k,ii)=0
2407                 do ik=i,j
2408                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2409                 enddo
2410                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2411                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2412              sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2413               else
2414                 if (odl(k,ii).le.dist_cut) then
2415                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2416                 else
2417 #ifdef OLDSIGMA
2418                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2419                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2420 #else
2421                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2422                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2423 #endif
2424                 endif
2425               endif
2426               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2427             else
2428               ii=ii+1
2429 !              l_homo(k,ii)=.false.
2430             endif
2431             enddo
2432           enddo
2433         lim_odl=ii
2434         endif
2435 !
2436 !     Theta, dihedral and SC retraints
2437 !
2438         if (waga_angle.gt.0.0d0) then
2439           do i = nnt_chain+3,nct_chain
2440             iii=i+chain_border1(1,igr)-1
2441             if (idomain(k,i).eq.0) then
2442 !               sigma_dih(k,i)=0.0
2443                cycle
2444             endif
2445             dih(k,iii)=phiref(i)
2446             sigma_dih(k,iii)= &
2447                (rescore(k,i)+rescore(k,i-1)+ &
2448                            rescore(k,i-2)+rescore(k,i-3))/4.0
2449 !            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2450 !     &       " sigma_dihed",sigma_dih(k,i)
2451             if (sigma_dih(k,iii).ne.0) &
2452              sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2453           enddo
2454 !          lim_dih=nct-nnt-2
2455         endif
2456
2457         if (waga_theta.gt.0.0d0) then
2458           do i = nnt_chain+2,nct_chain
2459              iii=i+chain_border1(1,igr)-1
2460              if (idomain(k,i).eq.0) then
2461 !              sigma_theta(k,i)=0.0
2462               cycle
2463              endif
2464              thetatpl(k,iii)=thetaref(i)
2465              sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
2466                               rescore(k,i-2))/3.0
2467              if (sigma_theta(k,iii).ne.0) &
2468              sigma_theta(k,iii)=1.0d0/ &
2469              (sigma_theta(k,iii)*sigma_theta(k,iii))
2470           enddo
2471         endif
2472
2473         if (waga_d.gt.0.0d0) then
2474           do i = nnt_chain,nct_chain
2475              iii=i+chain_border1(1,igr)-1
2476                if (itype(i,1).eq.10) cycle
2477                if (idomain(k,i).eq.0 ) then
2478 !                  sigma_d(k,i)=0.0
2479                   cycle
2480                endif
2481                xxtpl(k,iii)=xxref(i)
2482                yytpl(k,iii)=yyref(i)
2483                zztpl(k,iii)=zzref(i)
2484                sigma_d(k,iii)=rescore(k,i)
2485                if (sigma_d(k,iii).ne.0) &
2486                 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
2487 !               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
2488           enddo
2489         endif
2490       enddo ! l
2491       enddo ! ll
2492 !
2493 ! remove distance restraints not used in any model from the list
2494 ! shift data in all arrays
2495 !
2496 !      write (iout,*) "ii_old",ii_old
2497       if (waga_dist.ne.0.0d0) then
2498 #ifdef DEBUG
2499        write (iout,*) "Distance restraints from templates"
2500        do iii=1,lim_odl
2501        write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
2502         iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
2503         (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
2504         ki=1,constr_homology)
2505        enddo
2506 #endif
2507         ii=ii_old
2508         liiflag=.true.
2509         lfirst=.true.
2510         do i=nnt_chain,nct_chain-2
2511          do j=i+2,nct_chain
2512           ii=ii+1
2513 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2514 !     &            .and. distal.le.dist2_cut ) then
2515 !          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
2516 !          call flush(iout)
2517           if (ii_in_use(ii).eq.0.and.liiflag.or. &
2518           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2519             liiflag=.false.
2520             i10=ii
2521             if (lfirst) then
2522               lfirst=.false.
2523               iistart=ii
2524             else
2525               if(i10.eq.lim_odl) i10=i10+1
2526               do ki=0,i10-i01-1
2527                ires_homo(iistart+ki)=ires_homo(ki+i01)
2528                jres_homo(iistart+ki)=jres_homo(ki+i01)
2529                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2530                do k=1,constr_homology
2531                 odl(k,iistart+ki)=odl(k,ki+i01)
2532                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2533                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2534                enddo
2535               enddo
2536               iistart=iistart+i10-i01
2537             endif
2538           endif
2539           if (ii_in_use(ii).ne.0.and..not.liiflag) then
2540              i01=ii
2541              liiflag=.true.
2542           endif
2543          enddo
2544         enddo
2545         lim_odl=iistart-1
2546       endif
2547
2548       lll=lim_odl-ii_old
2549
2550       do i=2,nequiv(igr)
2551
2552         ichain=iequiv(i,igr)
2553
2554         do j=nnt_chain,nct_chain
2555           jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
2556           do k=1,constr_homology
2557             dih(k,jj)=dih(k,j)
2558             sigma_dih(k,jj)=sigma_dih(k,j)
2559             thetatpl(k,jj)=thetatpl(k,j)
2560             sigma_theta(k,jj)=sigma_theta(k,j)
2561             xxtpl(k,jj)=xxtpl(k,j)
2562             yytpl(k,jj)=yytpl(k,j)
2563             zztpl(k,jj)=zztpl(k,j)
2564             sigma_d(k,jj)=sigma_d(k,j)
2565           enddo
2566         enddo
2567
2568         jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
2569 !        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
2570         do j=ii_old+1,lim_odl
2571           ires_homo(j+lll)=ires_homo(j)+jj
2572           jres_homo(j+lll)=jres_homo(j)+jj
2573           do k=1,constr_homology
2574             odl(k,j+lll)=odl(k,j)
2575             sigma_odl(k,j+lll)=sigma_odl(k,j)
2576             l_homo(k,j+lll)=l_homo(k,j)
2577           enddo
2578         enddo
2579
2580         ii_old=ii_old+lll
2581         lim_odl=lim_odl+lll
2582
2583       enddo
2584
2585       ENDDO ! IGR
2586
2587       if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
2588       nres=nres_temp
2589       itype(:,1)=itype_temp(:)
2590
2591       return
2592    10 stop "Error in fragment file"
2593       end subroutine read_klapaucjusz
2594 !-----------------------------------------------------------------------------
2595       end module io