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