source/cluster/io_clust.F90

   1       module io_clust
   2 !-----------------------------------------------------------------------------
   3       use clust_data
   4       use io_units
   5       use geometry, only:int_from_cart,sc_loc_geom
   6       use io_config,only:readpdb,readpdb_template
   7 !      use names
   8       use io_base !, only: ilen
   9       use geometry_data, only: nres,c,nres_molec,boxxsize,boxysize,boxzsize
  10       use energy_data, only: nnt,nct,nss
  11       use control_data, only: lside
  12       implicit none
  13 !-----------------------------------------------------------------------------
  14 !
  15 !
  16 !-----------------------------------------------------------------------------
  17       contains
  18 !-----------------------------------------------------------------------------
  19 ! wrtclust.f
  20 !-----------------------------------------------------------------------------
  21       SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
  22
  23       use hc_, only: ioffset
  24       use control_data, only: lprint_cart,lprint_int,titel
  25       use geometry, only: int_from_cart1,nres
  26 !      implicit real*8 (a-h,o-z)
  27 !      include 'DIMENSIONS'
  28 !      include 'sizesclu.dat'
  29       integer,parameter :: num_in_line=5
  30       LOGICAL :: PRINTANG(max_cut)
  31       integer :: PRINTPDB(max_cut),printmol2(max_cut)
  32 !      include 'COMMON.CONTROL'
  33 !      include 'COMMON.HEADER'
  34 !      include 'COMMON.CHAIN'
  35 !      include 'COMMON.VAR'
  36 !      include 'COMMON.CLUSTER'
  37 !      include 'COMMON.IOUNITS'
  38 !      include 'COMMON.GEO'
  39 !      include 'COMMON.FREE'
  40 !      include 'COMMON.TEMPFAC'
  41       real(kind=8) :: rmsave(maxgr)
  42       CHARACTER(len=64) :: prefixp,NUMM,MUMM,EXTEN,extmol
  43       character(len=80) :: cfname
  44       character(len=8) :: ctemper
  45       DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,&
  46            MUMM /'000'/
  47 !      external ilen
  48       integer :: ncon,icut,ib
  49       integer :: i,ii,ii1,ii2,igr,ind1,ind2,ico,icon,&
  50                  irecord,nrecord,j,k,jj,ind,ncon_lim,ncon_out
  51       real(kind=8) :: temper,curr_dist,emin,qpart,boltz,&
  52                  ave_dim,amax_dim,emin1
  53
  54       if (allocated(tempfac)) deallocate(tempfac)
  55       allocate(tempfac(2,nres))
  56
  57       do i=1,64
  58         cfname(i:i)=" "
  59       enddo
  60 !      print *,"calling WRTCLUST",ncon
  61 !      write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
  62       rewind 80
  63       call flush(iout)
  64       temper=1.0d0/(beta_h(ib)*1.987d-3)
  65       if (temper.lt.100.0d0) then
  66         write(ctemper,'(f3.0)') temper
  67         ctemper(3:3)=" "
  68       else if (temper.lt.1000.0) then
  69         write (ctemper,'(f4.0)') temper
  70         ctemper(4:4)=" "
  71       else
  72         write (ctemper,'(f5.0)') temper
  73         ctemper(5:5)=" "
  74       endif
  75
  76       do i=1,ncon*(ncon-1)/2
  77         read (80) diss(i)
  78       enddo
  79       close(80,status='delete')
  80 !
  81 !  PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
  82 !
  83       ii1= index(intinname,'/')
  84       ii2=ii1
  85       ii1=ii1+1
  86       do while (ii2.gt.0)
  87         ii1=ii1+ii2
  88         ii2=index(intinname(ii1:),'/')
  89       enddo
  90       ii = ii1+index(intinname(ii1:),'.')-1
  91       if (ii.eq.0) then
  92         ii=ilen(intinname)
  93       else
  94         ii=ii-1
  95       endif
  96       prefixp=intinname(ii1:ii)
  97 !d    print *,icut,printang(icut),printpdb(icut),printmol2(icut)
  98 !d    print *,'ecut=',ecut
  99       WRITE (iout,100) NGR
 100       DO 19 IGR=1,NGR
 101       WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
 102       NRECORD=LICZ(IGR)/num_in_line
 103       IND1=1
 104       DO 63 IRECORD=1,NRECORD
 105       IND2=IND1+num_in_line-1
 106       WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
 107           totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
 108       IND1=IND2+1
 109    63 CONTINUE
 110       WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
 111          totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
 112       IND1=1
 113       ICON=list_conf(NCONF(IGR,1))
 114 !      WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
 115 ! 12/8/93 Estimation of "diameters" of the subsequent families.
 116       ave_dim=0.0
 117       amax_dim=0.0
 118 !      write (iout,*) "ecut",ecut
 119       do i=2,licz(igr)
 120         ii=nconf(igr,i)
 121         if (totfree(ii)-emin .gt. ecut) goto 10
 122         do j=1,i-1
 123           jj=nconf(igr,j)
 124           if (jj.eq.1) exit
 125           if (ii.lt.jj) then
 126             ind=ioffset(ncon,ii,jj)
 127           else
 128             ind=ioffset(ncon,jj,ii)
 129           endif
 130 !          write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
 131 !     &     " ind",ind
 132           call flush(iout)
 133           curr_dist=dabs(diss(ind)+0.0d0)
 134 !          write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
 135 !     &      list_conf(jj),curr_dist
 136           if (curr_dist .gt. amax_dim) amax_dim=curr_dist
 137           ave_dim=ave_dim+curr_dist**2
 138         enddo
 139       enddo
 140    10 if (licz(igr) .gt. 1) &
 141        ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
 142       write (iout,'(/A,F8.1,A,F8.1)') &
 143        'Max. distance in the family:',amax_dim,&
 144        '; average distance in the family:',ave_dim
 145       rmsave(igr)=0.0d0
 146       qpart=0.0d0
 147       do i=1,licz(igr)
 148         icon=nconf(igr,i)
 149         boltz=dexp(-totfree(icon))
 150         rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
 151         qpart=qpart+boltz
 152       enddo
 153       rmsave(igr)=rmsave(igr)/qpart
 154       write (iout,'(a,f5.2,a)') "Average RMSD",rmsave(igr)," A"
 155    19 CONTINUE
 156       WRITE (iout,400)
 157       WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
 158 !      print *,icut,printang(icut)
 159       IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
 160         emin=totfree_gr(1)
 161 !        print *,'emin',emin,' ngr',ngr
 162         if (lprint_cart) then
 163           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 164             //"K"//".x"
 165         else
 166           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 167             //"K"//".int"
 168         endif
 169         do igr=1,ngr
 170           icon=nconf(igr,1)
 171           if (totfree_gr(igr)-emin.le.ecut) then
 172             if (lprint_cart) then
 173               call cartout(igr,icon,totfree(icon)/beta_h(ib),&
 174                 totfree_gr(igr)/beta_h(ib),&
 175                 rmstb(icon),cfname)
 176             else
 177 !              print '(a)','calling briefout'
 178               do i=1,2*nres
 179                 do j=1,3
 180                   c(j,i)=allcart(j,i,icon)
 181                 enddo
 182               enddo
 183               call int_from_cart1(.false.)
 184 !el              call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),&
 185 !el                totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),&
 186 !el                jhpb_all(1,icon),cfname)
 187               call briefout(igr,totfree(icon)/beta_h(ib),&
 188                 totfree_gr(igr))
 189 !              print '(a)','exit briefout'
 190             endif
 191           endif
 192         enddo
 193         close(igeom)
 194       ENDIF
 195       IF (PRINTPDB(ICUT).gt.0) THEN
 196 ! Write out a number of conformations from each family in PDB format and
 197 ! create InsightII command file for their displaying in different colors
 198         cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 199           //"K_"//'ave'//exten
 200         write (iout,*) "cfname",cfname
 201         OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 202         write (ipdb,'(a,f8.2)') &
 203           "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
 204         close (ipdb)
 205         I=1
 206         ICON=NCONF(1,1)
 207         EMIN=totfree_gr(I)
 208         emin1=totfree(icon)
 209         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 210 !          write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
 211           write (NUMM,'(bz,i4.4)') i
 212           ncon_lim=min0(licz(i),printpdb(icut))
 213           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 214             //"K_"//numm(:ilen(numm))//exten
 215           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 216           write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5," AVE RMSD",0pf5.2)')&
 217            i,totfree_gr(i)/beta_h(ib),rmsave(i) !'
 218 ! Write conformations of the family i to PDB files
 219           ncon_out=1
 220           do while (ncon_out.lt.printpdb(icut) .and. &
 221            ncon_out.lt.licz(i).and. &
 222            totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 223             ncon_out=ncon_out+1
 224 !            write (iout,*) i,ncon_out,nconf(i,ncon_out),
 225 !     &        totfree(nconf(i,ncon_out)),emin1,ecut
 226           enddo
 227           write (iout,*) "ncon_out",ncon_out
 228           call flush(iout)
 229           do j=1,nres
 230             tempfac(1,j)=5.0d0
 231             tempfac(2,j)=5.0d0
 232           enddo
 233           do j=1,ncon_out
 234             icon=nconf(i,j)
 235             do ii=1,2*nres
 236               do k=1,3
 237                 c(k,ii)=allcart(k,ii,icon)
 238               enddo
 239             enddo
 240             call pdboutC(totfree(icon)/beta_h(ib),rmstb(icon),titel)
 241             write (ipdb,'("TER")')
 242           enddo
 243           close(ipdb)
 244 ! Average structures and structures closest to average
 245           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 246           //"K_"//'ave'//exten
 247           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',&
 248            position="APPEND")
 249           call ave_coord(i)
 250           write (ipdb,'(a,i5)') "REMARK CLUSTER",i
 251           call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 252           write (ipdb,'("TER")')
 253           call closest_coord(i)
 254           call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 255           write (ipdb,'("TER")')
 256           close (ipdb)
 257           I=I+1
 258           ICON=NCONF(I,1)
 259           emin1=totfree(icon)
 260         ENDDO
 261       ENDIF
 262       IF (printmol2(icut).gt.0) THEN
 263 ! Write out a number of conformations from each family in PDB format and
 264 ! create InsightII command file for their displaying in different colors
 265         I=1
 266         ICON=NCONF(1,1)
 267         EMIN=ENERGY(ICON)
 268         emin1=emin
 269         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 270           write (NUMM,'(bz,i4.4)') i
 271           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 272             //"K_"//numm(:ilen(numm))//extmol
 273           OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 274           ncon_out=1
 275           do while (ncon_out.lt.printmol2(icut) .and. &
 276            ncon_out.lt.licz(i).and. &
 277            totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 278             ncon_out=ncon_out+1
 279           enddo
 280           do j=1,ncon_out
 281             icon=nconf(i,j)
 282             do ii=1,2*nres
 283               do k=1,3
 284                 c(k,ii)=allcart(k,ii,icon)
 285               enddo
 286             enddo
 287             CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
 288           enddo
 289           CLOSE(imol2)
 290           I=I+1
 291           ICON=NCONF(I,1)
 292           emin1=totfree(icon)
 293         ENDDO
 294       ENDIF
 295   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
 296   200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,&
 297        ' CONTAINS ',I4,' CONFORMATION(S): ')
 298 ! 300 FORMAT ( 8(I4,F6.1))
 299   300 FORMAT (5(I4,1pe12.3))
 300   400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
 301   500 FORMAT (8(2I4,2X))
 302   600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
 303       RETURN
 304       END SUBROUTINE WRTCLUST
 305 !------------------------------------------------------------------------------
 306 !------------------------------------------------------------------------------
 307       subroutine ave_coord(igr)
 308
 309       use control_data, only:lside
 310       use regularize_, only:fitsq,matvec
 311 !      implicit none
 312 !      include 'DIMENSIONS'
 313 !      include 'sizesclu.dat'
 314 !      include 'COMMON.CONTROL'
 315 !      include 'COMMON.CLUSTER'
 316 !      include 'COMMON.CHAIN'
 317 !      include 'COMMON.INTERACT'
 318 !      include 'COMMON.VAR'
 319 !      include 'COMMON.TEMPFAC'
 320 !      include 'COMMON.IOUNITS'
 321       logical :: non_conv
 322       real(kind=8) :: przes(3),obrot(3,3)
 323       real(kind=8) :: xx(3,2*nres),yy(3,2*nres),csq(3,2*nres) !(3,maxres2)
 324       real(kind=8) :: eref
 325       integer :: i,ii,j,k,icon,jcon,igr
 326       real(kind=8) :: rms,boltz,qpart,cwork(3,2*nres),cref1(3,2*nres) !(3,maxres2)
 327 !      write (iout,*) "AVE_COORD: igr",igr
 328       jcon=nconf(igr,1)
 329       eref=totfree(jcon)
 330       boltz = dexp(-totfree(jcon)+eref)
 331       qpart=boltz
 332       do i=1,2*nres
 333         do j=1,3
 334           c(j,i)=allcart(j,i,jcon)*boltz
 335           cref1(j,i)=allcart(j,i,jcon)
 336           csq(j,i)=allcart(j,i,jcon)**2*boltz
 337         enddo
 338       enddo
 339       DO K=2,LICZ(IGR)
 340       jcon=nconf(igr,k)
 341       if (lside) then
 342         ii=0
 343         do i=nnt,nct
 344           ii=ii+1
 345           do j=1,3
 346             xx(j,ii)=allcart(j,i,jcon)
 347             yy(j,ii)=cref1(j,i)
 348           enddo
 349         enddo
 350         do i=nnt,nct
 351 !          if (itype(i).ne.10) then
 352             ii=ii+1
 353             do j=1,3
 354               xx(j,ii)=allcart(j,i+nres,jcon)
 355               yy(j,ii)=cref1(j,i+nres)
 356             enddo
 357 !          endif
 358         enddo
 359         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 360       else
 361         do i=nnt,nct
 362           do j=1,3
 363             cwork(j,i)=allcart(j,i,jcon)
 364           enddo
 365         enddo
 366         call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot &
 367              ,non_conv)
 368       endif
 369 !      write (iout,*) "rms",rms
 370 !      do i=1,3
 371 !        write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
 372 !      enddo
 373       if (rms.lt.0.0) then
 374         print *,'error, rms^2 = ',rms,icon,jcon
 375         stop
 376       endif
 377       if (non_conv) print *,non_conv,icon,jcon
 378       boltz=dexp(-totfree(jcon)+eref)
 379       qpart = qpart + boltz
 380       do i=1,2*nres
 381         do j=1,3
 382           xx(j,i)=allcart(j,i,jcon)
 383         enddo
 384       enddo
 385       call matvec(cwork,obrot,xx,2*nres)
 386       do i=1,2*nres
 387 !        write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
 388 !     &    (allcart(j,i,jcon),j=1,3)
 389         do j=1,3
 390           cwork(j,i)=cwork(j,i)+przes(j)
 391           c(j,i)=c(j,i)+cwork(j,i)*boltz
 392           csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
 393         enddo
 394       enddo
 395       ENDDO ! K
 396       do i=1,2*nres
 397         do j=1,3
 398           c(j,i)=c(j,i)/qpart
 399           csq(j,i)=csq(j,i)/qpart-c(j,i)**2
 400         enddo
 401 !        write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
 402       enddo
 403       do i=nnt,nct
 404         tempfac(1,i)=0.0d0
 405         tempfac(2,i)=0.0d0
 406         do j=1,3
 407           tempfac(1,i)=tempfac(1,i)+csq(j,i)
 408           tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
 409         enddo
 410         tempfac(1,i)=dsqrt(tempfac(1,i))
 411         tempfac(2,i)=dsqrt(tempfac(2,i))
 412       enddo
 413       return
 414       end subroutine ave_coord
 415 !------------------------------------------------------------------------------
 416       subroutine closest_coord(igr)
 417
 418       use regularize_, only:fitsq
 419 !      implicit none
 420 !      include 'DIMENSIONS'
 421 !      include 'sizesclu.dat'
 422 !      include 'COMMON.IOUNITS'
 423 !      include 'COMMON.CONTROL'
 424 !      include 'COMMON.CLUSTER'
 425 !      include 'COMMON.CHAIN'
 426 !      include 'COMMON.INTERACT'
 427 !      include 'COMMON.VAR'
 428       logical :: non_conv
 429       real(kind=8) :: przes(3),obrot(3,3)
 430       real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
 431       integer :: i,ii,j,k,icon,jcon,jconmin,igr
 432       real(kind=8) :: rms,rmsmin,cwork(3,2*nres)
 433       rmsmin=1.0d10
 434       jconmin=nconf(igr,1)
 435       DO K=1,LICZ(IGR)
 436       jcon=nconf(igr,k)
 437       if (lside) then
 438         ii=0
 439         do i=nnt,nct
 440           ii=ii+1
 441           do j=1,3
 442             xx(j,ii)=allcart(j,i,jcon)
 443             yy(j,ii)=c(j,i)
 444           enddo
 445         enddo
 446         do i=nnt,nct
 447 !          if (itype(i).ne.10) then
 448             ii=ii+1
 449             do j=1,3
 450               xx(j,ii)=allcart(j,i+nres,jcon)
 451               yy(j,ii)=c(j,i+nres)
 452             enddo
 453 !          endif
 454         enddo
 455         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 456       else
 457         do i=nnt,nct
 458           do j=1,3
 459             cwork(j,i)=allcart(j,i,jcon)
 460           enddo
 461         enddo
 462         call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot,&
 463              non_conv)
 464       endif
 465       if (rms.lt.0.0) then
 466         print *,'error, rms^2 = ',rms,icon,jcon
 467         stop
 468       endif
 469 !      write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
 470       if (non_conv) print *,non_conv,icon,jcon
 471       if (rms.lt.rmsmin) then
 472         rmsmin=rms
 473         jconmin=jcon
 474       endif
 475       ENDDO ! K
 476 !      write (iout,*) "rmsmin",rmsmin," rms",rms
 477       call flush(iout)
 478       do i=1,2*nres
 479         do j=1,3
 480           c(j,i)=allcart(j,i,jconmin)
 481         enddo
 482       enddo
 483       return
 484       end subroutine closest_coord
 485 !-----------------------------------------------------------------------------
 486 ! read_coords.F
 487 !-----------------------------------------------------------------------------
 488       subroutine read_coords(ncon,*)
 489
 490       use energy_data, only: ihpb,jhpb,max_ene
 491       use control_data, only: from_bx,from_cx
 492       use control, only: tcpu
 493 !      implicit none
 494 !      include "DIMENSIONS"
 495 !      include "sizesclu.dat"
 496 #ifdef MPI
 497       use MPI_data
 498       include "mpif.h"
 499       integer :: IERROR,ERRCODE !,STATUS(MPI_STATUS_SIZE)
 500 !      include "COMMON.MPI"
 501 #else
 502       use MPI_data, only: nprocs
 503 #endif
 504 !      include "COMMON.CONTROL"
 505 !      include "COMMON.CHAIN"
 506 !      include "COMMON.INTERACT"
 507 !      include "COMMON.IOUNITS"
 508 !      include "COMMON.VAR"
 509 !      include "COMMON.SBRIDGE"
 510 !      include "COMMON.GEO"
 511 !      include "COMMON.CLUSTER"
 512       character(len=3) :: liczba
 513       integer :: ncon
 514       integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,if,ib,&
 515         nn,nn1,inan
 516       integer :: ixdrf,iret,itmp
 517       real(kind=4) :: prec,reini,refree,rmsdev
 518       integer :: nrec,nlines,iscor,lenrec,lenrec_in
 519       real(kind=8) :: energ,t_acq !,tcpu
 520 !el      integer ilen,iroof
 521 !el      external ilen,iroof
 522       real(kind=8) :: rjunk
 523       integer :: ntot_all(0:nprocs-1) !(0:maxprocs-1)
 524       logical :: lerr
 525       real(kind=8) :: energia(0:max_ene),etot
 526       real(kind=4) :: csingle(3,2*nres+2)
 527       integer :: Previous,Next
 528       character(len=256) :: bprotfiles
 529 !      print *,"Processor",me," calls read_protein_data"
 530 #ifdef MPI
 531       if (me.eq.master) then
 532         Previous=MPI_PROC_NULL
 533       else
 534         Previous=me-1
 535       endif
 536       if (me.eq.nprocs-1) then
 537         Next=MPI_PROC_NULL
 538       else
 539         Next=me+1
 540       endif
 541 ! Set the scratchfile names
 542       write (liczba,'(bz,i3.3)') me
 543
 544       allocate(STATUS(MPI_STATUS_SIZE))
 545 #endif
 546 ! 1/27/05 AL Change stored coordinates to single precision and don't store
 547 !         energy components in the binary databases.
 548       lenrec=12*(nres+nct-nnt+1)+4*(2*nss+2)+16
 549       lenrec_in=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
 550 #ifdef DEBUG
 551       write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss", nss
 552       write (iout,*) "lenrec_in",lenrec_in
 553 #endif
 554       bprotfiles=scratchdir(:ilen(scratchdir))// &
 555              "/"//prefix(:ilen(prefix))//liczba//".xbin"
 556 ! EL
 557 ! allocate cluster arrays
 558       allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
 559       allocate(entfac(maxconf)) !(maxconf)
 560       allocate(rmstb(maxconf)) !(maxconf)
 561       allocate(allcart(3,2*nres,maxstr_proc)) !(3,maxres2,maxstr_proc)
 562       allocate(nss_all(maxstr_proc)) !(maxstr_proc)
 563       allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc)
 564       allocate(iscore(maxconf)) !(maxconf)
 565       nss_all=0
 566       ihpb_all=0
 567       jhpb_all=0
 568
 569 #ifdef CHUJ
 570       ICON=1
 571   123 continue
 572       if (from_cart .and. .not. from_bx .and. .not. from_cx) then
 573         if (efree) then
 574         read (intin,*,end=13,err=11) energy(icon),totfree(icon),&
 575           rmstb(icon),&
 576           nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
 577           i=1,nss_all(icon)),iscore(icon)
 578         else
 579         read (intin,*,end=13,err=11) energy(icon),rmstb(icon),&
 580           nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
 581           i=1,nss_all(icon)),iscore(icon)
 582         endif
 583         read (intin,'(8f10.5)',end=13,err=10) &
 584           ((allcart(j,i,icon),j=1,3),i=1,nres),&
 585           ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
 586         print *,icon,energy(icon),nss_all(icon),rmstb(icon)
 587       else
 588         read(intin,'(a80)',end=13,err=12) lineh
 589         read(lineh(:5),*,err=8) ic
 590         if (efree) then
 591         read(lineh(6:),*,err=8) energy(icon)
 592         else
 593         read(lineh(6:),*,err=8) energy(icon)
 594         endif
 595         goto 9
 596     8   ic=1
 597         print *,'error, assuming e=1d10',lineh
 598         energy(icon)=1d10
 599         nss=0
 600     9   continue
 601 !old        read(lineh(18:),*,end=13,err=11) nss_all(icon)
 602         ii = index(lineh(15:)," ")+15
 603         read(lineh(ii:),*,end=13,err=11) nss_all(icon)
 604         IF (NSS_all(icon).LT.9) THEN
 605           read (lineh(20:),*,end=102) &
 606           (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),&
 607           iscore(icon)
 608         ELSE
 609           read (lineh(20:),*,end=102) &
 610                  (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
 611           read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),&
 612             I=9,NSS_all(icon)),iscore(icon)
 613         ENDIF
 614
 615   102   continue
 616
 617         PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
 618         call read_angles(intin,*13)
 619         do i=1,nres
 620           phiall(i,icon)=phi(i)
 621           thetall(i,icon)=theta(i)
 622           alphall(i,icon)=alph(i)
 623           omall(i,icon)=omeg(i)
 624         enddo
 625       endif
 626       ICON=ICON+1
 627       GOTO 123
 628 !
 629 ! CALCULATE DISTANCES
 630 !
 631    10 print *,'something wrong with angles'
 632       goto 13
 633    11 print *,'something wrong with NSS',nss
 634       goto 13
 635    12 print *,'something wrong with header'
 636
 637    13 NCON=ICON-1
 638
 639 #endif
 640       call flush(iout)
 641       jj_old=1
 642       open (icbase,file=bprotfiles,status="unknown",&
 643          form="unformatted",access="direct",recl=lenrec)
 644 ! Read conformations from binary DA files (one per batch) and write them to
 645 ! a binary DA scratchfile.
 646       jj=0
 647       jjj=0
 648 #ifdef MPI
 649       write (liczba,'(bz,i3.3)') me
 650       IF (ME.EQ.MASTER) THEN
 651 ! Only the master reads the database; it'll send it to the other procs
 652 ! through a ring.
 653 #endif
 654         t_acq = tcpu()
 655         icount=0
 656
 657         if (from_bx) then
 658
 659           open (intin,file=intinname,status="old",form="unformatted",&
 660                   access="direct",recl=lenrec_in)
 661
 662         else if (from_cx) then
 663 #if (defined(AIX) && !defined(JUBL))
 664           call xdrfopen_(ixdrf,intinname, "r", iret)
 665 #else
 666           call xdrfopen(ixdrf,intinname, "r", iret)
 667 #endif
 668           prec=10000.0
 669           write (iout,*) "xdrfopen: iret",iret
 670           if (iret.eq.0) then
 671             write (iout,*) "Error: coordinate file ",&
 672              intinname(:ilen(intinname))," does not exist."
 673             call flush(iout)
 674 #ifdef MPI
 675             call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
 676 #endif
 677             stop
 678           endif
 679         else
 680           write (iout,*) "Error: coordinate format not specified"
 681           call flush(iout)
 682 #ifdef MPI
 683           call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
 684 #else
 685           stop
 686 #endif
 687         endif
 688
 689 !#define DEBUG
 690 #ifdef DEBUG
 691         write (iout,*) "Opening file ",intinname(:ilen(intinname))
 692         write (iout,*) "lenrec",lenrec_in
 693         call flush(iout)
 694 #endif
 695 !#undef DEBUG
 696 !        write (iout,*) "maxconf",maxconf
 697         i=0
 698         do while (.true.)
 699            i=i+1
 700 !el           if (i.gt.maxconf) then
 701 !el             write (iout,*) "Error: too many conformations ",&
 702 !el              "(",maxconf,") maximum."
 703 !#ifdef MPI
 704 !el             call MPI_Abort(MPI_COMM_WORLD,errcode,ierror)
 705 !#endif
 706 !el             stop
 707 !el           endif
 708 !          write (iout,*) "i",i
 709 !          call flush(iout)
 710           if (from_bx) then
 711             read(intin,err=101,end=101) &
 712              ((csingle(l,k),l=1,3),k=1,nres),&
 713              ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 714              nss,(ihpb(k),jhpb(k),k=1,nss),&
 715              energy(jj+1),&
 716              entfac(jj+1),rmstb(jj+1),iscor
 717              do j=1,2*nres
 718                do k=1,3
 719                  c(k,j)=csingle(k,j)
 720                enddo
 721              enddo
 722           else
 723 #if (defined(AIX) && !defined(JUBL))
 724             call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
 725             if (iret.eq.0) goto 101
 726             call xdrfint_(ixdrf, nss, iret)
 727             if (iret.eq.0) goto 101
 728             do j=1,nss
 729               call xdrfint_(ixdrf, ihpb(j), iret)
 730               if (iret.eq.0) goto 101
 731               call xdrfint_(ixdrf, jhpb(j), iret)
 732               if (iret.eq.0) goto 101
 733             enddo
 734             call xdrffloat_(ixdrf,reini,iret)
 735             if (iret.eq.0) goto 101
 736             call xdrffloat_(ixdrf,refree,iret)
 737             if (iret.eq.0) goto 101
 738             call xdrffloat_(ixdrf,rmsdev,iret)
 739             if (iret.eq.0) goto 101
 740             call xdrfint_(ixdrf,iscor,iret)
 741             if (iret.eq.0) goto 101
 742 #else
 743 !            write (iout,*) "calling xdrf3dfcoord"
 744             call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
 745 !            write (iout,*) "iret",iret
 746 !            call flush(iout)
 747             if (iret.eq.0) goto 101
 748             call xdrfint(ixdrf, nss, iret)
 749             write (iout,*) "iret",iret
 750 !            write (iout,*) "nss",nss,i,"TUTU"
 751             call flush(iout)
 752             if (iret.eq.0) goto 101
 753             do k=1,nss
 754               call xdrfint(ixdrf, ihpb(k), iret)
 755               if (iret.eq.0) goto 101
 756               call xdrfint(ixdrf, jhpb(k), iret)
 757               if (iret.eq.0) goto 101
 758 !            write(iout,*), ihpb(k),jhpb(k),"TUTU"
 759             enddo
 760             call xdrffloat(ixdrf,reini,iret)
 761             if (iret.eq.0) goto 101
 762 !            write(iout,*) "TUTU", reini
 763             call xdrffloat(ixdrf,refree,iret)
 764 !            write(iout,*) "TUTU", refree
 765             if (iret.eq.0) goto 101
 766             call xdrffloat(ixdrf,rmsdev,iret)
 767             if (iret.eq.0) goto 101
 768             call xdrfint(ixdrf,iscor,iret)
 769             if (iret.eq.0) goto 101
 770 #endif
 771             energy(jj+1)=reini
 772             entfac(jj+1)=refree
 773             rmstb(jj+1)=rmsdev
 774             do k=1,nres
 775               do l=1,3
 776                 c(l,k)=csingle(l,k)
 777               enddo
 778             enddo
 779             do k=nnt,nct
 780               do l=1,3
 781                 c(l,nres+k)=csingle(l,nres+k-nnt+1)
 782               enddo
 783             enddo
 784           endif
 785 #ifdef DEBUG
 786           write (iout,'(5hREAD ,i5,3f15.4,i10)') &
 787            jj+1,energy(jj+1),entfac(jj+1),&
 788            rmstb(jj+1),iscor
 789           write (iout,*) "Conformation",jjj+1,jj+1
 790           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 791           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 792           call flush(iout)
 793 #endif
 794           call add_new_cconf(jjj,jj,jj_old,icount,Next)
 795         enddo
 796   101   continue
 797         write (iout,*) i-1," conformations read from DA file ",&
 798           intinname(:ilen(intinname))
 799         write (iout,*) jj," conformations read so far"
 800         if (from_bx) then
 801           close(intin)
 802         else
 803 #if (defined(AIX) && !defined(JUBL))
 804           call xdrfclose_(ixdrf, iret)
 805 #else
 806           call xdrfclose(ixdrf, iret)
 807 #endif
 808         endif
 809 #ifdef MPI
 810 !#ifdef DEBUG
 811         write (iout,*) "jj_old",jj_old," jj",jj
 812 !#endif
 813         call write_and_send_cconf(icount,jj_old,jj,Next)
 814         call MPI_Send(0,1,MPI_INTEGER,Next,570,&
 815                    MPI_COMM_WORLD,IERROR)
 816         jj_old=jj+1
 817 #else
 818         call write_and_send_cconf(icount,jj_old,jj,Next)
 819 #endif
 820         t_acq = tcpu() - t_acq
 821 #ifdef MPI
 822         write (iout,*) "Processor",me,&
 823           " time for conformation read/send",t_acq
 824       ELSE
 825 ! A worker gets the confs from the master and sends them to its neighbor
 826         t_acq = tcpu()
 827         call receive_and_pass_cconf(icount,jj_old,jj,&
 828           Previous,Next)
 829         t_acq = tcpu() - t_acq
 830       ENDIF
 831 #endif
 832       ncon=jj
 833 !      close(icbase)
 834       close(intin)
 835
 836       write(iout,*)"A total of",ncon," conformations read."
 837
 838       allocate(enetb(1:max_ene,ncon)) !(max_ene,maxstr_proc)
 839 #ifdef MPI
 840 ! Check if everyone has the same number of conformations
 841       call MPI_Allgather(ncon,1,MPI_INTEGER,&
 842         ntot_all(0),1,MPI_INTEGER,MPI_Comm_World,IERROR)
 843       lerr=.false.
 844       do i=0,nprocs-1
 845         if (i.ne.me) then
 846             if (ncon.ne.ntot_all(i)) then
 847               write (iout,*) "Number of conformations at processor",i,&
 848                " differs from that at processor",me,&
 849                ncon,ntot_all(i)
 850               lerr = .true.
 851             endif
 852         endif
 853       enddo
 854       if (lerr) then
 855         write (iout,*)
 856         write (iout,*) "Number of conformations read by processors"
 857         write (iout,*)
 858         do i=0,nprocs-1
 859           write (iout,'(8i10)') i,ntot_all(i)
 860         enddo
 861         write (iout,*) "Calculation terminated."
 862         call flush(iout)
 863         return 1
 864       endif
 865       return
 866 #endif
 867  1111 write(iout,*) "Error opening coordinate file ",&
 868        intinname(:ilen(intinname))
 869       call flush(iout)
 870       return 1
 871       end subroutine read_coords
 872 !------------------------------------------------------------------------------
 873       subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
 874
 875       use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
 876       use energy_data, only: itel,itype,dsc,max_ene,molnum
 877       use control_data, only: symetr
 878       use geometry, only: int_from_cart1
 879 !      implicit none
 880 !      include "DIMENSIONS"
 881 !      include "sizesclu.dat"
 882 !      include "COMMON.CLUSTER"
 883 !      include "COMMON.CONTROL"
 884 !      include "COMMON.CHAIN"
 885 !      include "COMMON.INTERACT"
 886 !      include "COMMON.LOCAL"
 887 !      include "COMMON.IOUNITS"
 888 !      include "COMMON.NAMES"
 889 !      include "COMMON.VAR"
 890 !      include "COMMON.SBRIDGE"
 891 !      include "COMMON.GEO"
 892       integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
 893         nn,nn1,inan,Next,itj,chalen,mnum
 894       real(kind=8) :: etot,energia(0:max_ene)
 895       jjj=jjj+1
 896       chalen=int((nct-nnt+2)/symetr)
 897       call int_from_cart1(.false.)
 898       do j=nnt+1,nct-1
 899         mnum=molnum(j)
 900         write (iout,*) "Check atom",j,itype(j,mnum),vbld(j)
 901         if (mnum.ne.1) cycle
 902         if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
 903         if (itype(j-1,molnum(j-1)).eq.ntyp1_molec(molnum(j-1))) cycle
 904
 905         if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
 906          if (j.gt.2) then
 907           if (itel(j).ne.0 .and. itel(j-1).ne.0) then
 908           write (iout,*) "Conformation",jjj,jj+1,itype(j-1,molnum(j-1)),itype(j,mnum)
 909           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
 910        chalen
 911           write (iout,*) "The Cartesian geometry is:"
 912           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 913           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 914           write (iout,*) "The internal geometry is:"
 915           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 916           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 917           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 918           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 919           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 920           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 921           write (iout,*) &
 922             "This conformation WILL NOT be added to the database."
 923           return
 924           endif
 925          endif
 926         endif
 927       enddo
 928       do j=nnt,nct
 929         mnum=molnum(j)
 930         itj=itype(j,mnum)
 931         if (mnum.ne.1) cycle
 932         if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) &
 933             .and.  (vbld(nres+j)-dsc(iabs(itj))) &
 934                                   .gt.2.0d0) then
 935           write (iout,*) "Conformation",jjj,jj+1
 936           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
 937           write (iout,*) "The Cartesian geometry is:"
 938           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 939           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 940           write (iout,*) "The internal geometry is:"
 941           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 942           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 943           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 944           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 945           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 946           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 947           write (iout,*) &
 948             "This conformation WILL NOT be added to the database."
 949           return
 950         endif
 951       enddo
 952       do j=3,nres
 953         if (theta(j).le.0.0d0) then
 954           write (iout,*) &
 955             "Zero theta angle(s) in conformation",jjj,jj+1
 956           write (iout,*) "The Cartesian geometry is:"
 957           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 958           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 959           write (iout,*) "The internal geometry is:"
 960           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 961           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 962           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 963           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 964           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 965           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 966           write (iout,*) &
 967             "This conformation WILL NOT be added to the database."
 968           return
 969         endif
 970         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 971       enddo
 972       jj=jj+1
 973 #ifdef DEBUG
 974       write (iout,*) "Conformation",jjj,jj
 975       write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 976       write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 977       write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 978       write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 979       write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 980       write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 981       write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
 982       write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 983       write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 984       write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 985       write (iout,'(e15.5,16i5)') entfac(icount+1)
 986 !     &        iscore(icount+1,0)
 987 #endif
 988       icount=icount+1
 989       call store_cconf_from_file(jj,icount)
 990       if (icount.eq.maxstr_proc) then
 991 #ifdef DEBUG
 992         write (iout,* ) "jj_old",jj_old," jj",jj
 993 #endif
 994         call write_and_send_cconf(icount,jj_old,jj,Next)
 995         jj_old=jj+1
 996         icount=0
 997       endif
 998       return
 999       end subroutine add_new_cconf
1000 !------------------------------------------------------------------------------
1001       subroutine store_cconf_from_file(jj,icount)
1002
1003       use energy_data, only: ihpb,jhpb
1004 !      implicit none
1005 !      include "DIMENSIONS"
1006 !      include "sizesclu.dat"
1007 !      include "COMMON.CLUSTER"
1008 !      include "COMMON.CHAIN"
1009 !      include "COMMON.SBRIDGE"
1010 !      include "COMMON.INTERACT"
1011 !      include "COMMON.IOUNITS"
1012 !      include "COMMON.VAR"
1013       integer :: i,j,jj,icount
1014 ! Store the conformation that has been read in
1015       do i=1,2*nres
1016         do j=1,3
1017           allcart(j,i,icount)=c(j,i)
1018         enddo
1019       enddo
1020       nss_all(icount)=nss
1021       do i=1,nss
1022         ihpb_all(i,icount)=ihpb(i)
1023         jhpb_all(i,icount)=jhpb(i)
1024       enddo
1025       return
1026       end subroutine store_cconf_from_file
1027 !------------------------------------------------------------------------------
1028       subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1029
1030 !      implicit none
1031 !      include "DIMENSIONS"
1032 !      include "sizesclu.dat"
1033 #ifdef MPI
1034       use MPI_data
1035       include "mpif.h"
1036       integer :: IERROR
1037 !      include "COMMON.MPI"
1038 #endif
1039 !      include "COMMON.CHAIN"
1040 !      include "COMMON.SBRIDGE"
1041 !      include "COMMON.INTERACT"
1042 !      include "COMMON.IOUNITS"
1043 !      include "COMMON.CLUSTER"
1044 !      include "COMMON.VAR"
1045       integer :: icount,jj_old,jj,Next
1046 ! Write the structures to a scratch file
1047 #ifdef MPI
1048 ! Master sends the portion of conformations that have been read in to the neighbor
1049 #ifdef DEBUG
1050       write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1051       call flush(iout)
1052 #endif
1053       call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1054       call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1055           Next,571,MPI_COMM_WORLD,IERROR)
1056       call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1057           Next,572,MPI_COMM_WORLD,IERROR)
1058       call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1059           Next,573,MPI_COMM_WORLD,IERROR)
1060       call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1061           Next,577,MPI_COMM_WORLD,IERROR)
1062       call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1063           Next,579,MPI_COMM_WORLD,IERROR)
1064       call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1065           MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1066 #endif
1067       call dawrite_ccoords(jj_old,jj,icbase)
1068       return
1069       end subroutine write_and_send_cconf
1070 !------------------------------------------------------------------------------
1071 #ifdef MPI
1072       subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1073
1074       use MPI_data
1075 !      implicit none
1076 !      include "DIMENSIONS"
1077 !      include "sizesclu.dat"
1078       include "mpif.h"
1079       integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1080 !      include "COMMON.MPI"
1081 !      include "COMMON.CHAIN"
1082 !      include "COMMON.SBRIDGE"
1083 !      include "COMMON.INTERACT"
1084 !      include "COMMON.IOUNITS"
1085 !      include "COMMON.VAR"
1086 !      include "COMMON.GEO"
1087 !      include "COMMON.CLUSTER"
1088       integer :: i,j,k,icount,jj_old,jj,Previous,Next
1089       icount=1
1090 #ifdef DEBUG
1091       write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1092       call flush(iout)
1093 #endif
1094       do while (icount.gt.0)
1095       call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1096            STATUS,IERROR)
1097       call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1098            IERROR)
1099 #ifdef DEBUG
1100       write (iout,*) "Processor",me," icount",icount
1101 #endif
1102       if (icount.eq.0) return
1103       call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1104           Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1105       call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1106         Next,571,MPI_COMM_WORLD,IERROR)
1107       call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1108           Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1109       call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1110         Next,572,MPI_COMM_WORLD,IERROR)
1111       call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1112           Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1113       call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1114         Next,573,MPI_COMM_WORLD,IERROR)
1115       call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1116         Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1117       call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1118         Next,577,MPI_COMM_WORLD,IERROR)
1119       call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1120         Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1121       call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1122         Next,579,MPI_COMM_WORLD,IERROR)
1123       call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1124         MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1125       call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1126         MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1127       jj=jj_old+icount-1
1128       call dawrite_ccoords(jj_old,jj,icbase)
1129       jj_old=jj+1
1130 #ifdef DEBUG
1131       write (iout,*) "Processor",me," received",icount," conformations"
1132       do i=1,icount
1133         write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1134         write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1135         write (iout,'(e15.5,16i5)') entfac(i)
1136       enddo
1137 #endif
1138       enddo
1139       return
1140       end subroutine receive_and_pass_cconf
1141 #endif
1142 !------------------------------------------------------------------------------
1143       subroutine daread_ccoords(istart_conf,iend_conf)
1144
1145       use energy_data, only: dyn_ss
1146
1147 !      implicit none
1148 !      include "DIMENSIONS"
1149 !      include "sizesclu.dat"
1150 #ifdef MPI
1151       use MPI_data
1152       include "mpif.h"
1153 !      include "COMMON.MPI"
1154 #endif
1155 !      include "COMMON.CHAIN"
1156 !      include "COMMON.CLUSTER"
1157 !      include "COMMON.IOUNITS"
1158 !      include "COMMON.INTERACT"
1159 !      include "COMMON.VAR"
1160 !      include "COMMON.SBRIDGE"
1161 !      include "COMMON.GEO"
1162       integer :: istart_conf,iend_conf
1163       integer :: i,j,ij,ii,iii
1164       integer :: len
1165       character(len=16) :: form,acc
1166       character(len=32) :: nam
1167 !
1168 ! Read conformations off a DA scratchfile.
1169 !
1170 #ifdef DEBUG
1171       write (iout,*) "DAREAD_COORDS"
1172       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1173       inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1174       write (iout,*) "len=",len," form=",form," acc=",acc
1175       write (iout,*) "nam=",nam
1176       call flush(iout)
1177 #endif
1178       do ii=istart_conf,iend_conf
1179         ij = ii - istart_conf + 1
1180         iii=list_conf(ii)
1181 #ifdef DEBUG
1182         write (iout,*) "Reading binary file, record",iii," ii",ii
1183         call flush(iout)
1184 #endif
1185         if (dyn_ss) then
1186         read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
1187          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1188          entfac(ii),rmstb(ii)
1189         else
1190         read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1191          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1192          nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1193          entfac(ii),rmstb(ii)
1194         endif
1195
1196 #ifdef DEBUG
1197         write (iout,*) ii,iii,ij,entfac(ii)
1198         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1199         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1200           i=nnt+nres,nct+nres)
1201         write (iout,'(2e15.5)') entfac(ij)
1202         write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1203           jhpb_all(i,ij),i=1,nss)
1204         call flush(iout)
1205 #endif
1206       enddo
1207       return
1208       end subroutine daread_ccoords
1209 !------------------------------------------------------------------------------
1210       subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1211
1212 !      implicit none
1213 !      include "DIMENSIONS"
1214 !      include "sizesclu.dat"
1215       use energy_data, only: dyn_ss
1216
1217 #ifdef MPI
1218       use MPI_data
1219       include "mpif.h"
1220 !      include "COMMON.MPI"
1221 #endif
1222 !      include "COMMON.CHAIN"
1223 !      include "COMMON.INTERACT"
1224 !      include "COMMON.IOUNITS"
1225 !      include "COMMON.VAR"
1226 !      include "COMMON.SBRIDGE"
1227 !      include "COMMON.GEO"
1228 !      include "COMMON.CLUSTER"
1229       integer :: istart_conf,iend_conf
1230       integer :: i,j,ii,ij,iii,unit_out
1231       integer :: len
1232       character(len=16) :: form,acc
1233       character(len=32) :: nam
1234 !
1235 ! Write conformations to a DA scratchfile.
1236 !
1237 #ifdef DEBUG
1238       write (iout,*) "DAWRITE_COORDS"
1239       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1240       write (iout,*) "lenrec",lenrec
1241       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1242       write (iout,*) "len=",len," form=",form," acc=",acc
1243       write (iout,*) "nam=",nam
1244       call flush(iout)
1245 #endif
1246       do ii=istart_conf,iend_conf
1247         iii=list_conf(ii)
1248         ij = ii - istart_conf + 1
1249 #ifdef DEBUG
1250         write (iout,*) "Writing binary file, record",iii," ii",ii
1251         call flush(iout)
1252 #endif
1253        if (dyn_ss) then
1254         write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1255          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1256          entfac(ii),rmstb(ii)
1257         else
1258         write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
1259          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1260          nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1261          entfac(ii),rmstb(ii)
1262        endif
1263
1264 #ifdef DEBUG
1265         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1266         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1267          nct+nres)
1268         write (iout,'(2e15.5)') entfac(ij)
1269         write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1270          nss_all(ij))
1271         call flush(iout)
1272 #endif
1273       enddo
1274       return
1275       end subroutine dawrite_ccoords
1276 !-----------------------------------------------------------------------------
1277 ! readrtns.F
1278 !-----------------------------------------------------------------------------
1279       subroutine read_control
1280 !
1281 ! Read molecular data
1282 !
1283       use energy_data, only: rescale_mode,distchainmax,ipot,constr_homology,&
1284                  with_dihed_constr,read_homol_frag,out_template_coord,&
1285                  out_template_restr
1286       use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1287                  iscode,symetr,punch_dist,print_dist,nstart,nend,&
1288                  caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
1289                  r_cut_ele,nclust,tor_mode,scelemode,r_cut_mart,r_cut_ang,&
1290                  rlamb_mart,constr_dist
1291       use geometry_data, only:bordliptop,bordlipbot,&
1292           bufliptop,buflipbot,lipthick,lipbufthick
1293 !      implicit none
1294 !      include 'DIMENSIONS'
1295 !      include 'sizesclu.dat'
1296 !      include 'COMMON.IOUNITS'
1297 !      include 'COMMON.TIME1'
1298 !      include 'COMMON.SBRIDGE'
1299 !      include 'COMMON.CONTROL'
1300 !      include 'COMMON.CLUSTER'
1301 !      include 'COMMON.CHAIN'
1302 !      include 'COMMON.HEADER'
1303 !      include 'COMMON.FFIELD'
1304 !      include 'COMMON.FREE'
1305 !      include 'COMMON.INTERACT'
1306       character(len=320) :: controlcard !,ucase
1307 !#ifdef MPL
1308 !      include 'COMMON.INFO'
1309 !#endif
1310       integer :: i
1311
1312       read (INP,'(a80)') titel
1313       call card_concat(controlcard,.true.)
1314
1315       call readi(controlcard,'NRES',nres_molec(1),0)
1316       call readi(controlcard,'NRES_NUCL',nres_molec(2),0)
1317       call readi(controlcard,'NRES_CAT',nres_molec(5),0)
1318       nres=0
1319       do i=1,5
1320        nres=nres_molec(i)+nres
1321       enddo
1322       print *,"TU",nres_molec(:)
1323
1324 !      call alloc_clust_arrays
1325       allocate(rcutoff(max_cut+1)) !(max_cut+1)
1326       allocate(beta_h(maxT)) !(maxT)
1327       allocate(mult(nres)) !(maxres)
1328
1329
1330       call readi(controlcard,'RESCALE',rescale_mode,2)
1331       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1332       write (iout,*) "DISTCHAINMAX",distchainmax
1333       call readi(controlcard,'PDBOUT',outpdb,0)
1334       call readi(controlcard,'MOL2OUT',outmol2,0)
1335       refstr=(index(controlcard,'REFSTR').gt.0)
1336       write (iout,*) "REFSTR",refstr
1337       pdbref=(index(controlcard,'PDBREF').gt.0)
1338       iscode=index(controlcard,'ONE_LETTER')
1339       tree=(index(controlcard,'MAKE_TREE').gt.0)
1340       min_var=(index(controlcard,'MINVAR').gt.0)
1341       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1342       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1343       call readi(controlcard,'NCUT',ncut,1)
1344       call readi(controlcard,'SYM',symetr,1)
1345       write (iout,*) 'sym', symetr
1346       call readi(controlcard,'NSTART',nstart,0)
1347       call reada(controlcard,'BOXX',boxxsize,100.0d0)
1348       call reada(controlcard,'BOXY',boxysize,100.0d0)
1349       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1350       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
1351       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
1352       if (lipthick.gt.0.0d0) then
1353        bordliptop=(boxzsize+lipthick)/2.0
1354        bordlipbot=bordliptop-lipthick
1355       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
1356       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
1357       buflipbot=bordlipbot+lipbufthick
1358       bufliptop=bordliptop-lipbufthick
1359       if ((lipbufthick*2.0d0).gt.lipthick) &
1360        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
1361       endif !lipthick.gt.0
1362       write(iout,*) "bordliptop=",bordliptop
1363       write(iout,*) "bordlipbot=",bordlipbot
1364       write(iout,*) "bufliptop=",bufliptop
1365       write(iout,*) "buflipbot=",buflipbot
1366
1367       call readi(controlcard,'NCLUST',nclust,5)
1368 !      ions=index(controlcard,"IONS").gt.0
1369       call readi(controlcard,"SCELEMODE",scelemode,0)
1370       print *,"SCELE",scelemode
1371       call readi(controlcard,'TORMODE',tor_mode,0)
1372 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1373         write(iout,*) "torsional and valence angle mode",tor_mode
1374       call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
1375       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1376       write(iout,*) "R_CUT_ELE=",r_cut_ele
1377       call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
1378       call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
1379       call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
1380       call readi(controlcard,'NEND',nend,0)
1381       call reada(controlcard,'ECUT',ecut,10.0d0)
1382       call reada(controlcard,'PROB',prob_limit,0.99d0)
1383       write (iout,*) "Probability limit",prob_limit
1384       lgrp=(index(controlcard,'LGRP').gt.0)
1385       caonly=(index(controlcard,'CA_ONLY').gt.0)
1386       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1387       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1388       call readi(controlcard,'IOPT',iopt,2)
1389       lside = index(controlcard,"SIDE").gt.0
1390       efree = index(controlcard,"EFREE").gt.0
1391       call readi(controlcard,'NTEMP',nT,1)
1392       write (iout,*) "nT",nT
1393 !el      call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1394       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1395       write (iout,*) "nT",nT
1396       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1397       do i=1,nT
1398         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1399       enddo
1400       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1401       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1402       lprint_int=index(controlcard,"PRINT_INT") .gt.0
1403       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
1404       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
1405 !c      if (constr_homology) tole=dmax1(tole,1.5d0)
1406       write (iout,*) "with_homology_constr ",with_dihed_constr,&
1407       " CONSTR_HOMOLOGY",constr_homology
1408       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
1409       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
1410       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
1411       write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
1412       if (min_var) iopt=1
1413       return
1414       end subroutine read_control
1415 !-----------------------------------------------------------------------------
1416       subroutine read_constr_homology
1417       use MPI_data
1418       use energy_data
1419       use geometry_data
1420       use control_data
1421       use control, only:init_int_table,homology_partition
1422       use MD_data, only:iset
1423 !      implicit none
1424 !      include 'DIMENSIONS'
1425 !#ifdef MPI
1426 !      include 'mpif.h'
1427 !#endif
1428 !      include 'COMMON.SETUP'
1429 !      include 'COMMON.CONTROL'
1430 !      include 'COMMON.HOMOLOGY'
1431 !      include 'COMMON.CHAIN'
1432 !      include 'COMMON.IOUNITS'
1433 !      include 'COMMON.MD'
1434 !      include 'COMMON.QRESTR'
1435 !      include 'COMMON.GEO'
1436 !      include 'COMMON.INTERACT'
1437 !      include 'COMMON.NAMES'
1438 !      include 'COMMON.VAR'
1439 !
1440
1441 !     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1442 !    &                 dist_cut
1443 !     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1444 !    &    sigma_odl_temp(maxres,maxres,max_template)
1445       character*2 kic2
1446       character*24 model_ki_dist, model_ki_angle
1447       character*500 controlcard
1448       integer :: ki,i,ii,j,k,l
1449       integer, dimension (:), allocatable :: ii_in_use
1450       integer :: i_tmp,idomain_tmp,&
1451       irec,ik,iistart,nres_temp
1452 !      integer :: iset
1453 !      external :: ilen
1454       logical :: liiflag,lfirst
1455       integer :: i01,i10
1456 !
1457 !     FP - Nov. 2014 Temporary specifications for new vars
1458 !
1459       real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1460                        rescore3_tmp, dist_cut
1461       real(kind=8), dimension (:,:),allocatable :: rescore
1462       real(kind=8), dimension (:,:),allocatable :: rescore2
1463       real(kind=8), dimension (:,:),allocatable :: rescore3
1464       real(kind=8) :: distal
1465       character*24 tpl_k_rescore
1466       character*256 pdbfile
1467
1468 ! -----------------------------------------------------------------
1469 ! Reading multiple PDB ref structures and calculation of retraints
1470 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1471 ! FP (Nov., 2014)
1472 ! -----------------------------------------------------------------
1473 !
1474 !
1475 ! Alternative: reading from input
1476       call card_concat(controlcard,.true.)
1477       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1478       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1479       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1480       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1481       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1482       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1483       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1484       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1485       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1486       if(.not.read2sigma.and.start_from_model) then
1487           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1488            write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1489           start_from_model=.false.
1490           iranconf=(indpdb.le.0)
1491       endif
1492       if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1493          write(iout,*) 'START_FROM_MODELS is ON'
1494 !      if(start_from_model .and. rest) then
1495 !        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1496 !         write(iout,*) 'START_FROM_MODELS is OFF'
1497 !         write(iout,*) 'remove restart keyword from input'
1498 !        endif
1499 !      endif
1500       if (start_from_model) nmodel_start=constr_homology
1501       if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1502       if (homol_nset.gt.1)then
1503          call card_concat(controlcard,.true.)
1504          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1505          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1506 !          write(iout,*) "iset homology_weight "
1507           do i=1,homol_nset
1508            write(iout,*) i,waga_homology(i)
1509           enddo
1510          endif
1511          iset=mod(kolor,homol_nset)+1
1512       else
1513        iset=1
1514        waga_homology(1)=1.0
1515       endif
1516
1517 !d      write (iout,*) "nnt",nnt," nct",nct
1518 !d      call flush(iout)
1519
1520       if (read_homol_frag) then
1521         call read_klapaucjusz
1522       else
1523
1524       lim_odl=0
1525       lim_dih=0
1526 !
1527 !      write(iout,*) 'nnt=',nnt,'nct=',nct
1528 !
1529 !      do i = nnt,nct
1530 !        do k=1,constr_homology
1531 !         idomain(k,i)=0
1532 !        enddo
1533 !      enddo
1534        idomain=0
1535
1536 !      ii=0
1537 !      do i = nnt,nct-2
1538 !        do j=i+2,nct
1539 !        ii=ii+1
1540 !        ii_in_use(ii)=0
1541 !        enddo
1542 !      enddo
1543       ii_in_use=0
1544       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1545       if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
1546
1547       do k=1,constr_homology
1548
1549         read(inp,'(a)') pdbfile
1550         pdbfiles_chomo(k)=pdbfile
1551         if(me.eq.king .or. .not. out1file) &
1552          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1553         pdbfile(:ilen(pdbfile))
1554         open(ipdbin,file=pdbfile,status='old',err=33)
1555         goto 34
1556   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
1557         pdbfile(:ilen(pdbfile))
1558         stop
1559   34    continue
1560 !        print *,'Begin reading pdb data'
1561 !
1562 ! Files containing res sim or local scores (former containing sigmas)
1563 !
1564
1565         write(kic2,'(bz,i2.2)') k
1566
1567         tpl_k_rescore="template"//kic2//".sco"
1568         write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1569         unres_pdb=.false.
1570         nres_temp=nres
1571         write(iout,*) "read2sigma",read2sigma
1572
1573         if (read2sigma) then
1574           call readpdb_template(k)
1575         else
1576           call readpdb
1577         endif
1578         write(iout,*) "after readpdb"
1579         if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1580         nres_chomo(k)=nres
1581         nres=nres_temp
1582         if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1583         if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1584         if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1585         if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1586         if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1587         if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1588         if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1589         if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1590         if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1591         if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1592         if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1593         if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1594         if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1595         if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1596 !        if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1597         if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1598         if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1599         if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1600         if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1601 !        if(.not.allocated(distance)) allocate(distance(constr_homology))
1602 !        if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1603
1604
1605 !
1606 !     Distance restraints
1607 !
1608 !          ... --> odl(k,ii)
1609 ! Copy the coordinates from reference coordinates (?)
1610         do i=1,2*nres_chomo(k)
1611           do j=1,3
1612             c(j,i)=cref(j,i,1)
1613 !           write (iout,*) "c(",j,i,") =",c(j,i)
1614           enddo
1615         enddo
1616 !
1617 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1618 !
1619 !         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1620           open (ientin,file=tpl_k_rescore,status='old')
1621           if (nnt.gt.1) rescore(k,1)=0.0d0
1622           do irec=nnt,nct ! loop for reading res sim
1623             if (read2sigma) then
1624              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1625                                      rescore3_tmp,idomain_tmp
1626              i_tmp=i_tmp+nnt-1
1627              idomain(k,i_tmp)=idomain_tmp
1628              rescore(k,i_tmp)=rescore_tmp
1629              rescore2(k,i_tmp)=rescore2_tmp
1630              rescore3(k,i_tmp)=rescore3_tmp
1631              if (.not. out1file .or. me.eq.king)&
1632              write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
1633                            i_tmp,rescore2_tmp,rescore_tmp,&
1634                                      rescore3_tmp,idomain_tmp
1635             else
1636              idomain(k,irec)=1
1637              read (ientin,*,end=1401) rescore_tmp
1638
1639 !           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1640              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1641 !           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1642             endif
1643           enddo
1644  1401   continue
1645         close (ientin)
1646         if (waga_dist.ne.0.0d0) then
1647           ii=0
1648           do i = nnt,nct-2
1649             do j=i+2,nct
1650
1651               x12=c(1,i)-c(1,j)
1652               y12=c(2,i)-c(2,j)
1653               z12=c(3,i)-c(3,j)
1654               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1655 !              write (iout,*) k,i,j,distal,dist2_cut
1656
1657             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
1658                  .and. distal.le.dist2_cut ) then
1659
1660               ii=ii+1
1661               ii_in_use(ii)=1
1662               l_homo(k,ii)=.true.
1663
1664 !             write (iout,*) "k",k
1665 !             write (iout,*) "i",i," j",j," constr_homology",
1666 !    &                       constr_homology
1667               ires_homo(ii)=i
1668               jres_homo(ii)=j
1669               odl(k,ii)=distal
1670               if (read2sigma) then
1671                 sigma_odl(k,ii)=0
1672                 do ik=i,j
1673                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1674                 enddo
1675                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1676                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
1677               sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1678               else
1679                 if (odl(k,ii).le.dist_cut) then
1680                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1681                 else
1682 #ifdef OLDSIGMA
1683                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
1684                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1685 #else
1686                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
1687                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1688 #endif
1689                 endif
1690               endif
1691               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1692             else
1693 !              ii=ii+1
1694 !              l_homo(k,ii)=.false.
1695             endif
1696             enddo
1697           enddo
1698         lim_odl=ii
1699         endif
1700 !        write (iout,*) "Distance restraints set"
1701 !        call flush(iout)
1702 !
1703 !     Theta, dihedral and SC retraints
1704 !
1705         if (waga_angle.gt.0.0d0) then
1706 !         open (ientin,file=tpl_k_sigma_dih,status='old')
1707 !         do irec=1,maxres-3 ! loop for reading sigma_dih
1708 !            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1709 !            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1710 !            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1711 !    &                            sigma_dih(k,i+nnt-1)
1712 !         enddo
1713 !1402   continue
1714 !         close (ientin)
1715           do i = nnt+3,nct
1716             if (idomain(k,i).eq.0) then
1717                sigma_dih(k,i)=0.0
1718                cycle
1719             endif
1720             dih(k,i)=phiref(i) ! right?
1721 !           read (ientin,*) sigma_dih(k,i) ! original variant
1722 !             write (iout,*) "dih(",k,i,") =",dih(k,i)
1723 !             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1724 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1725 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1726 !    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1727
1728             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
1729                           rescore(k,i-2)+rescore(k,i-3))/4.0
1730 !            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1731 !           write (iout,*) "Raw sigmas for dihedral angle restraints"
1732 !           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1733 !           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1734 !                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1735 !   Instead of res sim other local measure of b/b str reliability possible
1736             if (sigma_dih(k,i).ne.0) &
1737             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1738 !           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1739           enddo
1740           lim_dih=nct-nnt-2
1741         endif
1742 !        write (iout,*) "Dihedral angle restraints set"
1743 !        call flush(iout)
1744
1745         if (waga_theta.gt.0.0d0) then
1746 !         open (ientin,file=tpl_k_sigma_theta,status='old')
1747 !         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1748 !            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1749 !            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1750 !    &                              sigma_theta(k,i+nnt-1)
1751 !         enddo
1752 !1403   continue
1753 !         close (ientin)
1754
1755           do i = nnt+2,nct ! right? without parallel.
1756 !         do i = i=1,nres ! alternative for bounds acc to readpdb?
1757 !         do i=ithet_start,ithet_end ! with FG parallel.
1758              if (idomain(k,i).eq.0) then
1759               sigma_theta(k,i)=0.0
1760               cycle
1761              endif
1762              thetatpl(k,i)=thetaref(i)
1763 !            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1764 !            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1765 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1766 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1767 !            read (ientin,*) sigma_theta(k,i) ! 1st variant
1768              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
1769                              rescore(k,i-2))/3.0
1770 !             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1771              if (sigma_theta(k,i).ne.0) &
1772              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1773
1774 !            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1775 !                             rescore(k,i-2) !  right expression ?
1776 !            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1777           enddo
1778         endif
1779 !        write (iout,*) "Angle restraints set"
1780 !        call flush(iout)
1781
1782         if (waga_d.gt.0.0d0) then
1783 !       open (ientin,file=tpl_k_sigma_d,status='old')
1784 !         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1785 !            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1786 !            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1787 !    &                          sigma_d(k,i+nnt-1)
1788 !         enddo
1789 !1404   continue
1790
1791           do i = nnt,nct ! right? without parallel.
1792 !         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1793 !         do i=loc_start,loc_end ! with FG parallel.
1794                if (itype(i,1).eq.10) cycle
1795                if (idomain(k,i).eq.0 ) then
1796                   sigma_d(k,i)=0.0
1797                   cycle
1798                endif
1799                xxtpl(k,i)=xxref(i)
1800                yytpl(k,i)=yyref(i)
1801                zztpl(k,i)=zzref(i)
1802 !              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1803 !              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1804 !              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1805 !              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1806                sigma_d(k,i)=rescore3(k,i) !  right expression ?
1807                if (sigma_d(k,i).ne.0) &
1808                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1809
1810 !              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1811 !              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1812 !              read (ientin,*) sigma_d(k,i) ! 1st variant
1813           enddo
1814         endif
1815       enddo
1816 !      write (iout,*) "SC restraints set"
1817 !      call flush(iout)
1818 !
1819 ! remove distance restraints not used in any model from the list
1820 ! shift data in all arrays
1821 !
1822 !      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
1823       if (waga_dist.ne.0.0d0) then
1824         ii=0
1825         liiflag=.true.
1826         lfirst=.true.
1827         do i=nnt,nct-2
1828          do j=i+2,nct
1829           ii=ii+1
1830 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1831 !     &            .and. distal.le.dist2_cut ) then
1832 !          write (iout,*) "i",i," j",j," ii",ii
1833 !          call flush(iout)
1834           if (ii_in_use(ii).eq.0.and.liiflag.or. &
1835           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
1836             liiflag=.false.
1837             i10=ii
1838             if (lfirst) then
1839               lfirst=.false.
1840               iistart=ii
1841             else
1842               if(i10.eq.lim_odl) i10=i10+1
1843               do ki=0,i10-i01-1
1844                ires_homo(iistart+ki)=ires_homo(ki+i01)
1845                jres_homo(iistart+ki)=jres_homo(ki+i01)
1846                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
1847                do k=1,constr_homology
1848                 odl(k,iistart+ki)=odl(k,ki+i01)
1849                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
1850                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
1851                enddo
1852               enddo
1853               iistart=iistart+i10-i01
1854             endif
1855           endif
1856           if (ii_in_use(ii).ne.0.and..not.liiflag) then
1857              i01=ii
1858              liiflag=.true.
1859           endif
1860          enddo
1861         enddo
1862         lim_odl=iistart-1
1863       endif
1864 !      write (iout,*) "Removing distances completed"
1865 !      call flush(iout)
1866       endif ! .not. klapaucjusz
1867
1868       if (constr_homology.gt.0) call homology_partition
1869       write (iout,*) "After homology_partition"
1870       call flush(iout)
1871       if (constr_homology.gt.0) call init_int_table
1872       write (iout,*) "After init_int_table"
1873       call flush(iout)
1874 !      endif ! .not. klapaucjusz
1875 !      endif
1876 !      if (constr_homology.gt.0) call homology_partition
1877 !      write (iout,*) "After homology_partition"
1878 !      call flush(iout)
1879 !      if (constr_homology.gt.0) call init_int_table
1880 !      write (iout,*) "After init_int_table"
1881 !      call flush(iout)
1882 !      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1883 !      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1884 !
1885 ! Print restraints
1886 !
1887       if (.not.out_template_restr) return
1888 !d      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1889       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1890        write (iout,*) "Distance restraints from templates"
1891        do ii=1,lim_odl
1892        write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
1893         ii,ires_homo(ii),jres_homo(ii),&
1894         (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
1895         ki=1,constr_homology)
1896        enddo
1897        write (iout,*) "Dihedral angle restraints from templates"
1898        do i=nnt+3,nct
1899         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
1900             (rad2deg*dih(ki,i),&
1901             rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1902        enddo
1903        write (iout,*) "Virtual-bond angle restraints from templates"
1904        do i=nnt+2,nct
1905         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
1906             (rad2deg*thetatpl(ki,i),&
1907             rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1908        enddo
1909        write (iout,*) "SC restraints from templates"
1910        do i=nnt,nct
1911         write(iout,'(i7,100(4f8.2,4x))') i,&
1912         (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
1913          1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1914        enddo
1915       endif
1916       return
1917       end subroutine read_constr_homology
1918 !----------------------------------------------------------------------------
1919       subroutine molread
1920 !
1921 ! Read molecular data.
1922 !
1923       use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc,&
1924                   deg2rad,phibound,crefjlee,cref,rad2deg
1925       use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1926 !                 wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1927 !                 wturn3,wturn4,wturn6,wvdwpp,weights
1928       use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1929                  indpdb,constr_dist,raw_psipred, with_theta_constr
1930       use geometry, only: chainbuild,alloc_geo_arrays
1931       use energy, only: alloc_ener_arrays
1932       use io_base, only: rescode
1933       use control, only: setup_var,init_int_table
1934       use conform_compar, only: contact
1935 !      implicit none
1936 !      include 'DIMENSIONS'
1937 !      include 'COMMON.IOUNITS'
1938 !      include 'COMMON.GEO'
1939 !      include 'COMMON.VAR'
1940 !      include 'COMMON.INTERACT'
1941 !      include 'COMMON.LOCAL'
1942 !      include 'COMMON.NAMES'
1943 !      include 'COMMON.CHAIN'
1944 !      include 'COMMON.FFIELD'
1945 !      include 'COMMON.SBRIDGE'
1946 !      include 'COMMON.HEADER'
1947 !      include 'COMMON.CONTROL'
1948 !      include 'COMMON.CONTACTS'
1949 !      include 'COMMON.TIME1'
1950 !#ifdef MPL
1951 !      include 'COMMON.INFO'
1952 !#endif
1953       character(len=4) ::  sequence(nres) !(maxres)
1954       character(len=800) :: weightcard,controlcard
1955 !      integer rescode
1956       real(kind=8) :: x(6*nres) !(maxvar)
1957       real(kind=8) :: ssscale
1958       real(kind=8) :: phihel,phibet,sigmahel,sigmabet,sumv,&
1959                       secprob(3,maxres)
1960       integer  :: itype_pdb(nres) !(maxres)
1961 !      logical seq_comp
1962       integer :: i,j,kkk,mnum
1963 !
1964 ! Body
1965 !
1966 !el      allocate(weights(n_ene))
1967       allocate(weights(max_ene))
1968       call alloc_geo_arrays
1969       call alloc_ener_arrays
1970 !-----------------------------
1971       allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1972       allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1973       allocate(itype(nres+2,5)) !(maxres)
1974       allocate(molnum(nres+1))
1975       allocate(itel(nres+2))
1976
1977       do i=1,2*nres+2
1978         do j=1,3
1979           c(j,i)=0
1980           dc(j,i)=0
1981         enddo
1982       enddo
1983       do j=1,5
1984       do i=1,nres+2
1985         itype(i,j)=0
1986         itel(i)=0
1987       enddo
1988       enddo
1989 !--------------------------
1990 ! Read weights of the subsequent energy terms.
1991       call card_concat(weightcard,.true.)
1992       call reada(weightcard,'WSC',wsc,1.0d0)
1993       call reada(weightcard,'WLONG',wsc,wsc)
1994       call reada(weightcard,'WSCP',wscp,1.0d0)
1995       call reada(weightcard,'WELEC',welec,1.0D0)
1996       call reada(weightcard,'WVDWPP',wvdwpp,welec)
1997       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1998       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1999       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
2000       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
2001       call reada(weightcard,'WTURN3',wturn3,1.0D0)
2002       call reada(weightcard,'WTURN4',wturn4,1.0D0)
2003       call reada(weightcard,'WTURN6',wturn6,1.0D0)
2004       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
2005       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
2006       call reada(weightcard,'WBOND',wbond,1.0D0)
2007       call reada(weightcard,'WTOR',wtor,1.0D0)
2008       call reada(weightcard,'WTORD',wtor_d,1.0D0)
2009       call reada(weightcard,'WANG',wang,1.0D0)
2010       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
2011       call reada(weightcard,'SCAL14',scal14,0.4D0)
2012       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
2013       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
2014       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
2015       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
2016       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
2017       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
2018       if (index(weightcard,'SOFT').gt.0) ipot=6
2019       call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
2020       call reada(weightcard,'WELPP',welpp,0.0d0)
2021       call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
2022       call reada(weightcard,'WELPSB',welpsb,0.0D0)
2023       call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
2024       call reada(weightcard,'WELSB',welsb,0.0D0)
2025       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
2026       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
2027       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
2028       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
2029 !      print *,"KUR..",wtor_nucl
2030       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
2031       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
2032       call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
2033 ! 12/1/95 Added weight for the multi-body term WCORR
2034       call reada(weightcard,'WCORRH',wcorr,1.0D0)
2035       call reada(weightcard,'WSCBASE',wscbase,0.0D0)
2036       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
2037       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
2038       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
2039       call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
2040       call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
2041       call reada(weightcard,"D0CM",d0cm,3.78d0)
2042       call reada(weightcard,"AKCM",akcm,15.1d0)
2043       call reada(weightcard,"AKTH",akth,11.0d0)
2044       call reada(weightcard,"AKCT",akct,12.0d0)
2045       call reada(weightcard,"V1SS",v1ss,-1.08d0)
2046       call reada(weightcard,"V2SS",v2ss,7.61d0)
2047       call reada(weightcard,"V3SS",v3ss,13.7d0)
2048       call reada(weightcard,"EBR",ebr,-5.50D0)
2049       call reada(weightcard,"ATRISS",atriss,0.301D0)
2050       call reada(weightcard,"BTRISS",btriss,0.021D0)
2051       call reada(weightcard,"CTRISS",ctriss,1.001D0)
2052       call reada(weightcard,"DTRISS",dtriss,1.001D0)
2053       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
2054       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
2055       call reada(weightcard,"LIPSCALE",lipscale,1.0D0)
2056       call reada(weightcard,"WTL",wliptran,1.0D0)
2057
2058       call reada(weightcard,"HT",Ht,0.0D0)
2059       if (dyn_ss) then
2060        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
2061         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
2062         akcm=akcm/wsc*ssscale
2063         akth=akth/wsc*ssscale
2064         akct=akct/wsc*ssscale
2065         v1ss=v1ss/wsc*ssscale
2066         v2ss=v2ss/wsc*ssscale
2067         v3ss=v3ss/wsc*ssscale
2068       else
2069         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
2070       endif
2071
2072       if (wcorr4.gt.0.0d0) wcorr=wcorr4
2073       weights(1)=wsc
2074       weights(2)=wscp
2075       weights(3)=welec
2076       weights(4)=wcorr
2077       weights(5)=wcorr5
2078       weights(6)=wcorr6
2079       weights(7)=wel_loc
2080       weights(8)=wturn3
2081       weights(9)=wturn4
2082       weights(10)=wturn6
2083       weights(11)=wang
2084       weights(12)=wscloc
2085       weights(13)=wtor
2086       weights(14)=wtor_d
2087       weights(15)=wstrain
2088       weights(16)=wvdwpp
2089       weights(17)=wbond
2090       weights(18)=scal14
2091 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
2092       weights(21)=wsccor
2093           weights(26)=wvdwpp_nucl
2094           weights(27)=welpp
2095           weights(28)=wvdwpsb
2096           weights(29)=welpsb
2097           weights(30)=wvdwsb
2098           weights(31)=welsb
2099           weights(32)=wbond_nucl
2100           weights(33)=wang_nucl
2101           weights(34)=wsbloc
2102           weights(35)=wtor_nucl
2103           weights(36)=wtor_d_nucl
2104           weights(37)=wcorr_nucl
2105           weights(38)=wcorr3_nucl
2106           weights(41)=wcatcat
2107           weights(42)=wcatprot
2108           weights(46)=wscbase
2109           weights(47)=wscpho
2110           weights(48)=wpeppho
2111           weights(49)=wpeppho
2112           weights(50)=wcatnucl
2113           weights(56)=wcat_tran
2114           weights(22)=wliptran
2115
2116       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
2117         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
2118         wturn4,wturn6,wsccor
2119    10 format (/'Energy-term weights (unscaled):'// &
2120        'WSCC=   ',f10.6,' (SC-SC)'/ &
2121        'WSCP=   ',f10.6,' (SC-p)'/ &
2122        'WELEC=  ',f10.6,' (p-p electr)'/ &
2123        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
2124        'WBOND=  ',f10.6,' (stretching)'/ &
2125        'WANG=   ',f10.6,' (bending)'/ &
2126        'WSCLOC= ',f10.6,' (SC local)'/ &
2127        'WTOR=   ',f10.6,' (torsional)'/ &
2128        'WTORD=  ',f10.6,' (double torsional)'/ &
2129        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
2130        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
2131        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
2132        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
2133        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
2134        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
2135        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
2136        'WTURN6= ',f10.6,' (turns, 6th order)'/ &
2137        'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
2138
2139       if (wcorr4.gt.0.0d0) then
2140         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
2141          'between contact pairs of peptide groups'
2142         write (iout,'(2(a,f5.3/))') &
2143         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
2144         'Range of quenching the correlation terms:',2*delt_corr
2145       else if (wcorr.gt.0.0d0) then
2146         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
2147          'between contact pairs of peptide groups'
2148       endif
2149       write (iout,'(a,f8.3)') &
2150         'Scaling factor of 1,4 SC-p interactions:',scal14
2151       write (iout,'(a,f8.3)') &
2152         'General scaling factor of SC-p interactions:',scalscp
2153       r0_corr=cutoff_corr-delt_corr
2154       do i=1,20
2155         aad(i,1)=scalscp*aad(i,1)
2156         aad(i,2)=scalscp*aad(i,2)
2157         bad(i,1)=scalscp*bad(i,1)
2158         bad(i,2)=scalscp*bad(i,2)
2159       enddo
2160
2161       call flush(iout)
2162       print *,'indpdb=',indpdb,' pdbref=',pdbref
2163
2164 ! Read sequence if not taken from the pdb file.
2165       if (iscode.gt.0) then
2166         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
2167       else
2168         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
2169       endif
2170       print *,nres_molec(:),nres
2171 ! Convert sequence to numeric code
2172       do i=1,nres_molec(1)
2173         mnum=1
2174         molnum(i)=1
2175         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2176       enddo
2177       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
2178         mnum=2
2179         molnum(i)=2
2180         write (iout,*),i,sequence(i)
2181         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2182
2183       enddo
2184       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
2185         mnum=5
2186         molnum(i)=5
2187         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2188       enddo
2189       print *,nres
2190       print '(20i4)',(itype(i,molnum(i)),i=1,nres)
2191
2192       do i=1,nres-1
2193 #ifdef PROCOR
2194         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
2195             .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
2196 #else
2197         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
2198 #endif
2199           itel(i)=0
2200 #ifdef PROCOR
2201         else if (iabs(itype(i+1,1)).ne.20) then
2202 #else
2203         else if (iabs(itype(i,1)).ne.20) then
2204 #endif
2205           itel(i)=1
2206         else
2207           itel(i)=2
2208         endif
2209       enddo
2210       write (iout,*) "ITEL"
2211       do i=1,nres-1
2212         write (iout,*) i,itype(i,molnum(i)),itel(i)
2213       enddo
2214
2215       print *,'Call Read_Bridge.'
2216       call read_bridge
2217       nnt=1
2218       nct=nres
2219       print *,'NNT=',NNT,' NCT=',NCT
2220       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
2221       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
2222       if (nstart.lt.nnt) nstart=nnt
2223       if (nend.gt.nct .or. nend.eq.0) nend=nct
2224       write (iout,*) "nstart",nstart," nend",nend
2225       nres0=nres
2226 !      if (pdbref) then
2227 !        read(inp,'(a)') pdbfile
2228 !        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
2229 !        open(ipdbin,file=pdbfile,status='old',err=33)
2230 !        goto 34
2231 !  33    write (iout,'(a)') 'Error opening PDB file.'
2232 !        stop
2233 !  34    continue
2234 !        print *,'Begin reading pdb data'
2235 !        call readpdb
2236 !        print *,'Finished reading pdb data'
2237 !        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
2238 !        do i=1,nres
2239 !          itype_pdb(i)=itype(i)
2240 !        enddo
2241 !        close (ipdbin)
2242 !        write (iout,'(a,i3)') 'nsup=',nsup
2243 !        nstart_seq=nnt
2244 !        if (nsup.le.(nct-nnt+1)) then
2245 !          do i=0,nct-nnt+1-nsup
2246 !            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
2247 !              nstart_seq=nnt+i
2248 !              goto 111
2249 !            endif
2250 !          enddo
2251 !          write (iout,'(a)')
2252 !     &            'Error - sequences to be superposed do not match.'
2253 !          stop
2254 !        else
2255 !          do i=0,nsup-(nct-nnt+1)
2256 !            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
2257 !     &      then
2258 !              nstart_sup=nstart_sup+i
2259 !              nsup=nct-nnt+1
2260 !              goto 111
2261 !            endif
2262 !          enddo
2263 !          write (iout,'(a)')
2264 !     &            'Error - sequences to be superposed do not match.'
2265 !        endif
2266 !  111   continue
2267 !        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
2268 !     &                 ' nstart_seq=',nstart_seq
2269 !      endif
2270             if (with_dihed_constr) then
2271
2272       read (inp,*) ndih_constr
2273       if (ndih_constr.gt.0) then
2274         raw_psipred=.false.
2275 !C        read (inp,*) ftors
2276 !C        write (iout,*) 'FTORS',ftors
2277 !C ftors is the force constant for torsional quartic constrains
2278         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),&
2279                      i=1,ndih_constr)
2280         write (iout,*)&
2281         'There are',ndih_constr,' constraints on phi angles.'
2282         do i=1,ndih_constr
2283           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),&
2284        ftors(i)
2285         enddo
2286         do i=1,ndih_constr
2287           phi0(i)=deg2rad*phi0(i)
2288           drange(i)=deg2rad*drange(i)
2289         enddo
2290       else if (ndih_constr.lt.0) then
2291         raw_psipred=.true.
2292         call card_concat(controlcard,.true.)
2293         call reada(controlcard,"PHIHEL",phihel,50.0D0)
2294         call reada(controlcard,"PHIBET",phibet,180.0D0)
2295         call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
2296         call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
2297         call reada(controlcard,"WDIHC",wdihc,0.591d0)
2298         write (iout,*) "Weight of the dihedral restraint term",wdihc
2299         read(inp,'(9x,3f7.3)')&
2300           (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
2301         write (iout,*) "The secprob array"
2302         do i=nnt,nct
2303           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
2304         enddo
2305         ndih_constr=0
2306         do i=nnt+3,nct
2307           if (itype(i-3,molnum(i-3)).ne.ntyp1 .and. itype(i-2,molnum(i-2)).ne.ntyp1&
2308          .and. itype(i-1,molnum(i-1)).ne.ntyp1 .and. itype(i,molnum(i)).ne.ntyp1) then
2309             ndih_constr=ndih_constr+1
2310             idih_constr(ndih_constr)=i
2311             sumv=0.0d0
2312             do j=1,3
2313               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
2314               sumv=sumv+vpsipred(j,ndih_constr)
2315             enddo
2316             do j=1,3
2317               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
2318             enddo
2319             phibound(1,ndih_constr)=phihel*deg2rad
2320             phibound(2,ndih_constr)=phibet*deg2rad
2321             sdihed(1,ndih_constr)=sigmahel*deg2rad
2322             sdihed(2,ndih_constr)=sigmabet*deg2rad
2323           endif
2324         enddo
2325         write (iout,*)&
2326         'There are',ndih_constr,&
2327         ' bimodal restraints on gamma angles.'
2328         do i=1,ndih_constr
2329           write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,&
2330            restyp(itype(idih_constr(i)-2,molnum(idih_constr(i)-2)),&
2331                  molnum(idih_constr(i)-2)),idih_constr(i)-2,&
2332            restyp(itype(idih_constr(i)-1,molnum(idih_constr(i)-1)),&
2333                  molnum(idih_constr(i)-2)),idih_constr(i)-1,&
2334            phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,&
2335            phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,&
2336            (vpsipred(j,i),j=1,3)
2337         enddo
2338
2339       endif ! endif ndif_constr.gt.0
2340       endif ! with_dihed_constr
2341       if (with_theta_constr) then
2342 !C with_theta_constr is keyword allowing for occurance of theta constrains
2343       read (inp,*) ntheta_constr
2344 !C ntheta_constr is the number of theta constrains
2345       if (ntheta_constr.gt.0) then
2346 !C        read (inp,*) ftors
2347         read (inp,*) (itheta_constr(i),theta_constr0(i),&
2348        theta_drange(i),for_thet_constr(i),&
2349        i=1,ntheta_constr)
2350 !C the above code reads from 1 to ntheta_constr
2351 !C itheta_constr(i) residue i for which is theta_constr
2352 !C theta_constr0 the global minimum value
2353 !C theta_drange is range for which there is no energy penalty
2354 !C for_thet_constr is the force constant for quartic energy penalty
2355 !C E=k*x**4
2356          write (iout,*)&
2357         'There are',ntheta_constr,' constraints on phi angles.'
2358          do i=1,ntheta_constr
2359           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),&
2360          theta_drange(i),&
2361          for_thet_constr(i)
2362          enddo
2363 !C        endif
2364         do i=1,ntheta_constr
2365          theta_constr0(i)=deg2rad*theta_constr0(i)
2366           theta_drange(i)=deg2rad*theta_drange(i)
2367         enddo
2368 !C        if(me.eq.king.or..not.out1file)
2369 !C     &   write (iout,*) 'FTORS',ftors
2370 !C        do i=1,ntheta_constr
2371 !C          ii = itheta_constr(i)
2372 !C          thetabound(1,ii) = phi0(i)-drange(i)
2373 !C          thetabound(2,ii) = phi0(i)+drange(i)
2374 !C        enddo
2375       endif ! ntheta_constr.gt.0
2376       endif! with_theta_constr
2377       if (constr_homology.gt.0) then
2378 !c        write (iout,*) "About to call read_constr_homology"
2379 !c        call flush(iout)
2380         call read_constr_homology
2381 !c        write (iout,*) "Exit read_constr_homology"
2382 !c        call flush(iout)
2383         if (indpdb.gt.0 .or. pdbref) then
2384           kkk=1
2385           do i=1,2*nres
2386             do j=1,3
2387               c(j,i)=crefjlee(j,i)
2388               cref(j,i,kkk)=crefjlee(j,i)
2389             enddo
2390           enddo
2391         endif
2392 #ifdef DEBUG
2393         write (iout,*) "Array C"
2394         do i=1,nres
2395           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
2396            (c(j,i+nres),j=1,3)
2397         enddo
2398         write (iout,*) "Array Cref"
2399         do i=1,nres
2400           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),&
2401            (cref(j,i+nres),j=1,3)
2402         enddo
2403 #endif
2404 #ifdef DEBUG
2405        call int_from_cart1(.false.)
2406        call sc_loc_geom(.false.)
2407        do i=1,nres
2408          thetaref(i)=theta(i)
2409          phiref(i)=phi(i)
2410          write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
2411        enddo
2412        do i=1,nres-1
2413          do j=1,3
2414            dc(j,i)=c(j,i+1)-c(j,i)
2415            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2416          enddo
2417        enddo
2418        do i=2,nres-1
2419          do j=1,3
2420            dc(j,i+nres)=c(j,i+nres)-c(j,i)
2421            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2422          enddo
2423        enddo
2424 #endif
2425       else
2426         homol_nset=0
2427       endif
2428
2429
2430       write(iout,*)"przed ini_int_tab"
2431       call init_int_table
2432       write(iout,*)"po ini_int_tab"
2433       write(iout,*)"przed setup var"
2434       call setup_var
2435       write(iout,*)"po setup var"
2436       if (ns.gt.0.and.dyn_ss) then
2437           do i=nss+1,nhpb
2438             ihpb(i-nss)=ihpb(i)
2439             jhpb(i-nss)=jhpb(i)
2440             forcon(i-nss)=forcon(i)
2441             dhpb(i-nss)=dhpb(i)
2442           enddo
2443           nhpb=nhpb-nss
2444           nss=0
2445           link_start=1
2446           link_end=nhpb
2447 !          call hpb_partition
2448           do i=1,ns
2449             dyn_ss_mask(iss(i))=.true.
2450           enddo
2451       endif
2452
2453       write (iout,*) "molread: REFSTR",refstr
2454       if (refstr) then
2455         if (.not.pdbref) then
2456           call read_angles(inp,*38)
2457           goto 39
2458    38     write (iout,'(a)') 'Error reading reference structure.'
2459 #ifdef MPL
2460           call mp_stopall(Error_Msg)
2461 #else
2462           stop 'Error reading reference structure'
2463 #endif
2464    39     call chainbuild
2465           nstart_sup=nnt
2466           nstart_seq=nnt
2467           nsup=nct-nnt+1
2468           kkk=1
2469           do i=1,2*nres
2470             do j=1,3
2471               cref(j,i,kkk)=c(j,i)
2472             enddo
2473           enddo
2474         endif
2475         call contact(.true.,ncont_ref,icont_ref)
2476       endif
2477       return
2478       end subroutine molread
2479 !-----------------------------------------------------------------------------
2480       subroutine openunits
2481 !      implicit none
2482 !      include 'DIMENSIONS'
2483       use control_data, only: from_cx,from_bx,from_cart
2484 #ifdef MPI
2485       use MPI_data
2486       include "mpif.h"
2487       character(len=3) :: liczba
2488 !      include "COMMON.MPI"
2489 #endif
2490 !      include 'COMMON.IOUNITS'
2491 !      include 'COMMON.CONTROL'
2492       integer :: lenpre,lenpot !,ilen
2493 !      external ilen
2494       character(len=16) :: cformat,cprint
2495 !      character(len=16) ucase
2496       integer :: lenint,lenout
2497       call getenv('INPUT',prefix)
2498       call getenv('OUTPUT',prefout)
2499       call getenv('INTIN',prefintin)
2500       call getenv('COORD',cformat)
2501       call getenv('PRINTCOOR',cprint)
2502       call getenv('SCRATCHDIR',scratchdir)
2503       from_bx=.true.
2504       from_cx=.false.
2505       if (index(ucase(cformat),'CX').gt.0) then
2506         from_cx=.true.
2507         from_bx=.false.
2508       endif
2509       from_cart=.true.
2510       lenpre=ilen(prefix)
2511       lenout=ilen(prefout)
2512       lenint=ilen(prefintin)
2513 ! Get the names and open the input files
2514       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
2515 #ifdef MPI
2516       write (liczba,'(bz,i3.3)') me
2517       outname=prefout(:lenout)//'_clust.out_'//liczba
2518 #else
2519       outname=prefout(:lenout)//'_clust.out'
2520 #endif
2521       if (from_bx) then
2522         intinname=prefintin(:lenint)//'.bx'
2523       else if (from_cx) then
2524         intinname=prefintin(:lenint)//'.cx'
2525       else
2526         intinname=prefintin(:lenint)//'.int'
2527       endif
2528       rmsname=prefintin(:lenint)//'.rms'
2529       open (jplot,file=prefout(:ilen(prefout))//'.tex',&
2530              status='unknown')
2531       open (jrms,file=rmsname,status='unknown')
2532       open(iout,file=outname,status='unknown')
2533 ! Get parameter filenames and open the parameter files.
2534       call getenv('BONDPAR',bondname)
2535       open (ibond,file=bondname,status='old')
2536       call getenv('THETPAR',thetname)
2537       open (ithep,file=thetname,status='old')
2538       call getenv('ROTPAR',rotname)
2539       open (irotam,file=rotname,status='old')
2540       call getenv('TORPAR',torname)
2541       open (itorp,file=torname,status='old')
2542       call getenv('TORDPAR',tordname)
2543       open (itordp,file=tordname,status='old')
2544       call getenv('FOURIER',fouriername)
2545       open (ifourier,file=fouriername,status='old')
2546       call getenv('ELEPAR',elename)
2547       open (ielep,file=elename,status='old')
2548       call getenv('SIDEPAR',sidename)
2549       open (isidep,file=sidename,status='old')
2550       call getenv('SIDEP',sidepname)
2551       open (isidep1,file=sidepname,status="old")
2552       call getenv('SCCORPAR',sccorname)
2553       open (isccor,file=sccorname,status="old")
2554       call getenv('BONDPAR_NUCL',bondname_nucl)
2555       open (ibond_nucl,file=bondname_nucl,status='old')
2556       call getenv('THETPAR_NUCL',thetname_nucl)
2557       open (ithep_nucl,file=thetname_nucl,status='old')
2558       call getenv('ROTPAR_NUCL',rotname_nucl)
2559       open (irotam_nucl,file=rotname_nucl,status='old')
2560       call getenv('TORPAR_NUCL',torname_nucl)
2561       open (itorp_nucl,file=torname_nucl,status='old')
2562       call getenv('TORDPAR_NUCL',tordname_nucl)
2563       open (itordp_nucl,file=tordname_nucl,status='old')
2564       call getenv('SIDEPAR_NUCL',sidename_nucl)
2565       open (isidep_nucl,file=sidename_nucl,status='old')
2566       call getenv('SIDEPAR_SCBASE',sidename_scbase)
2567       open (isidep_scbase,file=sidename_scbase,status='old')
2568       call getenv('PEPPAR_PEPBASE',pepname_pepbase)
2569       open (isidep_pepbase,file=pepname_pepbase,status='old')
2570       call getenv('SCPAR_PHOSPH',pepname_scpho)
2571       open (isidep_scpho,file=pepname_scpho,status='old')
2572       call getenv('PEPPAR_PHOSPH',pepname_peppho)
2573       open (isidep_peppho,file=pepname_peppho,status='old')
2574
2575
2576       call getenv('LIPTRANPAR',liptranname)
2577       open (iliptranpar,file=liptranname,status='old')
2578       call getenv('TUBEPAR',tubename)
2579       open (itube,file=tubename,status='old')
2580       call getenv('IONPAR',ionname)
2581       open (iion,file=ionname,status='old')
2582       call getenv('IONPAR_NUCL',ionnuclname)
2583       open (iionnucl,file=ionnuclname,status='old')
2584       call getenv('IONPAR_TRAN',iontranname)
2585       open (iiontran,file=iontranname,status='old')
2586
2587
2588 #ifndef OLDSCP
2589 !
2590 ! 8/9/01 In the newest version SCp interaction constants are read from a file
2591 ! Use -DOLDSCP to use hard-coded constants instead.
2592 !
2593       call getenv('SCPPAR',scpname)
2594       open (iscpp,file=scpname,status='old')
2595       call getenv('SCPPAR_NUCL',scpname_nucl)
2596       open (iscpp_nucl,file=scpname_nucl,status='old')
2597
2598 #endif
2599       return
2600       end subroutine openunits
2601 !-----------------------------------------------------------------------------
2602 ! geomout.F
2603 !-----------------------------------------------------------------------------
2604       subroutine pdboutC(etot,rmsd,tytul)
2605
2606       use energy_data, only: ihpb,jhpb,itype,molnum
2607       use energy, only:boxshift,to_box
2608 !      implicit real*8 (a-h,o-z)
2609 !      include 'DIMENSIONS'
2610 !      include 'COMMON.CONTROL'
2611 !      include 'COMMON.CHAIN'
2612 !      include 'COMMON.INTERACT'
2613 !      include 'COMMON.NAMES'
2614 !      include 'COMMON.IOUNITS'
2615 !      include 'COMMON.HEADER'
2616 !      include 'COMMON.SBRIDGE'
2617 !      include 'COMMON.TEMPFAC'
2618       character(len=50) :: tytul
2619       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
2620                                                'G','H','I','J'/)
2621       integer :: ica(nres),k,iti1,ki
2622       real(kind=8) :: etot,rmsd
2623       integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
2624       real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
2625       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
2626         ' ENERGY ',etot,' RMS ',rmsd
2627       iatom=0
2628       ichain=1
2629       ires=0
2630 !      boxxxshift(1)=int(c(1,nnt)/boxxsize)
2631 !      if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
2632 !      boxxxshift(2)=int(c(2,nnt)/boxzsize)
2633 !      if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
2634 !      boxxxshift(3)=int(c(3,nnt)/boxzsize)
2635 !      if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
2636 !      do i=1,3
2637 !        if (boxxxshift(i).gt.2) boxxxshift(i)=2
2638 !        if (boxxxshift(i).lt.-2) boxxxshift(i)=-2
2639 !
2640 !      enddo
2641       do j=1,3
2642       cbeg(j)=c(j,nnt)
2643       enddo
2644       call to_box(cbeg(1),cbeg(2),cbeg(3))
2645       boxxxx(1)=boxxsize
2646       boxxxx(2)=boxysize
2647       boxxxx(3)=boxzsize
2648
2649       do i=nnt,nct
2650         mnum=molnum(i)
2651         iti=itype(i,mnum)
2652         if (i.ne.nct) then
2653         iti1=itype(i+1,mnum)
2654         mnum1=molnum(i+1)
2655         if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle
2656         endif
2657         if (iti.eq.ntyp1_molec(mnum)) then
2658           ichain=ichain+1
2659           ires=0
2660           write (ipdb,'(a)') 'TER'
2661         else
2662         ires=ires+1
2663 !         if (molnum(i).ge.1) then
2664           if (i.le.3) then
2665            ki=2
2666           else
2667           if (itype(i,mnum).ne.ntyp1_molec(mnum)) then
2668           ki=ki+1
2669           endif
2670          endif
2671 !        endif
2672 !          print *,"tu2",i,k
2673           if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
2674           if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
2675         do j=1,3
2676         cbeg(j)=c(j,ki)
2677         enddo
2678         iatom=iatom+1
2679         ica(i)=iatom
2680         call to_box(c(1,i),c(2,i),c(3,i))
2681
2682         DO k=1,3
2683         Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k))
2684         c(k,i)=cbeg(k)+Rdist(k)
2685         ENDDO
2686       if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then
2687         call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres))
2688         DO k=1,3
2689         Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k))
2690         c(k,i+nres)=cbeg(k)+Rdist(k)
2691         ENDDO
2692        endif
2693         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
2694            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
2695         if ((iti.ne.10).and.(mnum.ne.5)) then
2696           iatom=iatom+1
2697           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
2698             ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
2699         endif
2700         endif
2701       enddo
2702       write (ipdb,'(a)') 'TER'
2703       do i=nnt,nct-1
2704         mnum=molnum(i)
2705         mnum1=molnum(i+1)
2706         if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
2707         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2708           write (ipdb,30) ica(i),ica(i+1)
2709         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2710           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2711         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
2712           write (ipdb,30) ica(i),ica(i)+1
2713         endif
2714       enddo
2715       if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
2716         write (ipdb,30) ica(nct),ica(nct)+1
2717       endif
2718       do i=1,nss
2719         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2720       enddo
2721       write (ipdb,'(a6)') 'ENDMDL'
2722   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2723   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2724   30  FORMAT ('CONECT',8I5)
2725       return
2726       end subroutine pdboutC
2727 !-----------------------------------------------------------------------------
2728       subroutine cartout(igr,i,etot,free,rmsd,plik)
2729 !     implicit real*8 (a-h,o-z)
2730 !     include 'DIMENSIONS'
2731 !     include 'sizesclu.dat'
2732 !     include 'COMMON.IOUNITS'
2733 !     include 'COMMON.CHAIN'
2734 !     include 'COMMON.VAR'
2735 !     include 'COMMON.LOCAL'
2736 !     include 'COMMON.INTERACT'
2737 !     include 'COMMON.NAMES'
2738 !     include 'COMMON.GEO'
2739 !     include 'COMMON.CLUSTER'
2740       integer :: igr,i,j,k
2741       real(kind=8) :: etot,free,rmsd
2742       character(len=80) :: plik
2743       open (igeom,file=plik,position='append')
2744       write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
2745       write (igeom,'(i4,$)') &
2746         nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
2747       write (igeom,'(i10)') iscore(i)
2748       write (igeom,'(8f10.5)') &
2749         ((allcart(k,j,i),k=1,3),j=1,nres),&
2750         ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
2751       return
2752       end subroutine cartout
2753 !------------------------------------------------------------------------------
2754       subroutine read_klapaucjusz
2755       use energy_data
2756       use control_data, only:unres_pdb
2757       use geometry_data, only:theta,thetaref,phi,phiref,&
2758       xxref,yyref,zzref
2759       implicit none
2760 !     include 'DIMENSIONS'
2761 !#ifdef MPI
2762 !     include 'mpif.h'
2763 !#endif
2764 !     include 'COMMON.SETUP'
2765 !     include 'COMMON.CONTROL'
2766 !     include 'COMMON.HOMOLOGY'
2767 !     include 'COMMON.CHAIN'
2768 !     include 'COMMON.IOUNITS'
2769 !     include 'COMMON.MD'
2770 !     include 'COMMON.GEO'
2771 !     include 'COMMON.INTERACT'
2772 !     include 'COMMON.NAMES'
2773       character(len=256) fragfile
2774       integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2775                          ii_in_use
2776       integer, dimension(:,:), allocatable :: iresclust,inclust
2777       integer :: nclust
2778
2779       character(len=2) :: kic2
2780       character(len=24) :: model_ki_dist, model_ki_angle
2781       character(len=500) :: controlcard
2782       integer :: ki, i, j, jj,k, l, i_tmp,&
2783       idomain_tmp,&
2784       ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2785       i01,i10,nnt_chain,nct_chain
2786       real(kind=8) :: distal
2787       logical :: lprn = .true.
2788       integer :: nres_temp
2789 !      integer :: ilen
2790 !      external :: ilen
2791       logical :: liiflag,lfirst
2792
2793       real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2794       real(kind=8), dimension (:,:), allocatable  :: rescore
2795       real(kind=8), dimension (:,:), allocatable :: rescore2
2796       character(len=24) :: tpl_k_rescore
2797       character(len=256) :: pdbfile
2798
2799 !
2800 ! For new homol impl
2801 !
2802 !     include 'COMMON.VAR'
2803 !
2804 !      write (iout,*) "READ_KLAPAUCJUSZ"
2805 !      print *,"READ_KLAPAUCJUSZ"
2806 !      call flush(iout)
2807       call getenv("FRAGFILE",fragfile)
2808       write (iout,*) "Opening", fragfile
2809       call flush(iout)
2810       open(ientin,file=fragfile,status="old",err=10)
2811 !      write (iout,*) " opened"
2812 !      call flush(iout)
2813
2814       sigma_theta=0.0
2815       sigma_d=0.0
2816       sigma_dih=0.0
2817       l_homo = .false.
2818
2819       nres_temp=nres
2820       itype_temp(:)=itype(:,1)
2821       ii=0
2822       lim_odl=0
2823
2824 !      write (iout,*) "Entering loop"
2825 !      call flush(iout)
2826
2827       DO IGR = 1,NCHAIN_GROUP
2828
2829 !      write (iout,*) "igr",igr
2830       call flush(iout)
2831       read(ientin,*) constr_homology,nclust
2832       if (start_from_model) then
2833         nmodel_start=constr_homology
2834       else
2835         nmodel_start=0
2836       endif
2837
2838       ii_old=lim_odl
2839
2840       ichain=iequiv(1,igr)
2841       nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2842       nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2843 !      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2844
2845 ! Read pdb files
2846       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2847       if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2848       do k=1,constr_homology
2849         read(ientin,'(a)') pdbfile
2850         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2851         pdbfile(:ilen(pdbfile))
2852         pdbfiles_chomo(k)=pdbfile
2853         open(ipdbin,file=pdbfile,status='old',err=33)
2854         goto 34
2855   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
2856         pdbfile(:ilen(pdbfile))
2857         stop
2858   34    continue
2859         unres_pdb=.false.
2860 !        nres_temp=nres
2861         call readpdb_template(k)
2862         nres_chomo(k)=nres
2863 !        nres=nres_temp
2864         do i=1,nres
2865           rescore(k,i)=0.2d0
2866           rescore2(k,i)=1.0d0
2867         enddo
2868       enddo
2869 ! Read clusters
2870       do i=1,nclust
2871         read(ientin,*) ninclust(i),nresclust(i)
2872         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2873         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2874       enddo
2875 !
2876 ! Loop over clusters
2877 !
2878       do l=1,nclust
2879         do ll = 1,ninclust(l)
2880
2881         k = inclust(ll,l)
2882 !        write (iout,*) "l",l," ll",ll," k",k
2883         do i=1,nres
2884           idomain(k,i)=0
2885         enddo
2886         do i=1,nresclust(l)
2887           if (nnt.gt.1)  then
2888             idomain(k,iresclust(i,l)+1) = 1
2889           else
2890             idomain(k,iresclust(i,l)) = 1
2891           endif
2892         enddo
2893 !
2894 !     Distance restraints
2895 !
2896 !          ... --> odl(k,ii)
2897 ! Copy the coordinates from reference coordinates (?)
2898 !        nres_temp=nres
2899         nres=nres_chomo(k)
2900         do i=1,2*nres
2901           do j=1,3
2902             c(j,i)=chomo(j,i,k)
2903 !           write (iout,*) "c(",j,i,") =",c(j,i)
2904           enddo
2905         enddo
2906         call int_from_cart(.true.,.false.)
2907         call sc_loc_geom(.false.)
2908         do i=1,nres
2909           thetaref(i)=theta(i)
2910           phiref(i)=phi(i)
2911         enddo
2912 !        nres=nres_temp
2913         if (waga_dist.ne.0.0d0) then
2914           ii=ii_old
2915 !          do i = nnt,nct-2
2916           do i = nnt_chain,nct_chain-2
2917 !            do j=i+2,nct
2918             do j=i+2,nct_chain
2919
2920               x12=c(1,i)-c(1,j)
2921               y12=c(2,i)-c(2,j)
2922               z12=c(3,i)-c(3,j)
2923               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2924 !              write (iout,*) k,i,j,distal,dist2_cut
2925
2926             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2927                  .and. distal.le.dist2_cut ) then
2928
2929               ii=ii+1
2930               ii_in_use(ii)=1
2931               l_homo(k,ii)=.true.
2932
2933 !             write (iout,*) "k",k
2934 !             write (iout,*) "i",i," j",j," constr_homology",
2935 !     &                       constr_homology
2936               ires_homo(ii)=i+chain_border1(1,igr)-1
2937               jres_homo(ii)=j+chain_border1(1,igr)-1
2938               odl(k,ii)=distal
2939               if (read2sigma) then
2940                 sigma_odl(k,ii)=0
2941                 do ik=i,j
2942                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2943                 enddo
2944                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2945                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2946              sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2947               else
2948                 if (odl(k,ii).le.dist_cut) then
2949                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2950                 else
2951 #ifdef OLDSIGMA
2952                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2953                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2954 #else
2955                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2956                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2957 #endif
2958                 endif
2959               endif
2960               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2961             else
2962               ii=ii+1
2963 !              l_homo(k,ii)=.false.
2964             endif
2965             enddo
2966           enddo
2967         lim_odl=ii
2968         endif
2969 !
2970 !     Theta, dihedral and SC retraints
2971 !
2972         if (waga_angle.gt.0.0d0) then
2973           do i = nnt_chain+3,nct_chain
2974             iii=i+chain_border1(1,igr)-1
2975             if (idomain(k,i).eq.0) then
2976 !               sigma_dih(k,i)=0.0
2977                cycle
2978             endif
2979             dih(k,iii)=phiref(i)
2980             sigma_dih(k,iii)= &
2981                (rescore(k,i)+rescore(k,i-1)+ &
2982                            rescore(k,i-2)+rescore(k,i-3))/4.0
2983 !            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2984 !     &       " sigma_dihed",sigma_dih(k,i)
2985             if (sigma_dih(k,iii).ne.0) &
2986              sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2987           enddo
2988 !          lim_dih=nct-nnt-2
2989         endif
2990
2991         if (waga_theta.gt.0.0d0) then
2992           do i = nnt_chain+2,nct_chain
2993              iii=i+chain_border1(1,igr)-1
2994              if (idomain(k,i).eq.0) then
2995 !              sigma_theta(k,i)=0.0
2996               cycle
2997              endif
2998              thetatpl(k,iii)=thetaref(i)
2999              sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
3000                               rescore(k,i-2))/3.0
3001              if (sigma_theta(k,iii).ne.0) &
3002              sigma_theta(k,iii)=1.0d0/ &
3003              (sigma_theta(k,iii)*sigma_theta(k,iii))
3004           enddo
3005         endif
3006
3007         if (waga_d.gt.0.0d0) then
3008           do i = nnt_chain,nct_chain
3009              iii=i+chain_border1(1,igr)-1
3010                if (itype(i,1).eq.10) cycle
3011                if (idomain(k,i).eq.0 ) then
3012 !                  sigma_d(k,i)=0.0
3013                   cycle
3014                endif
3015                xxtpl(k,iii)=xxref(i)
3016                yytpl(k,iii)=yyref(i)
3017                zztpl(k,iii)=zzref(i)
3018                sigma_d(k,iii)=rescore(k,i)
3019                if (sigma_d(k,iii).ne.0) &
3020                 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
3021 !               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3022           enddo
3023         endif
3024       enddo ! l
3025       enddo ! ll
3026 !
3027 ! remove distance restraints not used in any model from the list
3028 ! shift data in all arrays
3029 !
3030 !      write (iout,*) "ii_old",ii_old
3031       if (waga_dist.ne.0.0d0) then
3032 #ifdef DEBUG
3033        write (iout,*) "Distance restraints from templates"
3034        do iii=1,lim_odl
3035        write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
3036         iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
3037         (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
3038         ki=1,constr_homology)
3039        enddo
3040 #endif
3041         ii=ii_old
3042         liiflag=.true.
3043         lfirst=.true.
3044         do i=nnt_chain,nct_chain-2
3045          do j=i+2,nct_chain
3046           ii=ii+1
3047 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3048 !     &            .and. distal.le.dist2_cut ) then
3049 !          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
3050 !          call flush(iout)
3051           if (ii_in_use(ii).eq.0.and.liiflag.or. &
3052           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3053             liiflag=.false.
3054             i10=ii
3055             if (lfirst) then
3056               lfirst=.false.
3057               iistart=ii
3058             else
3059               if(i10.eq.lim_odl) i10=i10+1
3060               do ki=0,i10-i01-1
3061                ires_homo(iistart+ki)=ires_homo(ki+i01)
3062                jres_homo(iistart+ki)=jres_homo(ki+i01)
3063                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3064                do k=1,constr_homology
3065                 odl(k,iistart+ki)=odl(k,ki+i01)
3066                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3067                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3068                enddo
3069               enddo
3070               iistart=iistart+i10-i01
3071             endif
3072           endif
3073           if (ii_in_use(ii).ne.0.and..not.liiflag) then
3074              i01=ii
3075              liiflag=.true.
3076           endif
3077          enddo
3078         enddo
3079         lim_odl=iistart-1
3080       endif
3081
3082       lll=lim_odl-ii_old
3083
3084       do i=2,nequiv(igr)
3085
3086         ichain=iequiv(i,igr)
3087
3088         do j=nnt_chain,nct_chain
3089           jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
3090           do k=1,constr_homology
3091             dih(k,jj)=dih(k,j)
3092             sigma_dih(k,jj)=sigma_dih(k,j)
3093             thetatpl(k,jj)=thetatpl(k,j)
3094             sigma_theta(k,jj)=sigma_theta(k,j)
3095             xxtpl(k,jj)=xxtpl(k,j)
3096             yytpl(k,jj)=yytpl(k,j)
3097             zztpl(k,jj)=zztpl(k,j)
3098             sigma_d(k,jj)=sigma_d(k,j)
3099           enddo
3100         enddo
3101
3102         jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
3103 !        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
3104         do j=ii_old+1,lim_odl
3105           ires_homo(j+lll)=ires_homo(j)+jj
3106           jres_homo(j+lll)=jres_homo(j)+jj
3107           do k=1,constr_homology
3108             odl(k,j+lll)=odl(k,j)
3109             sigma_odl(k,j+lll)=sigma_odl(k,j)
3110             l_homo(k,j+lll)=l_homo(k,j)
3111           enddo
3112         enddo
3113
3114         ii_old=ii_old+lll
3115         lim_odl=lim_odl+lll
3116
3117       enddo
3118
3119       ENDDO ! IGR
3120
3121       if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
3122       nres=nres_temp
3123       itype(:,1)=itype_temp(:)
3124
3125       return
3126    10 stop "Error in fragment file"
3127       end subroutine read_klapaucjusz
3128
3129 !-----------------------------------------------------------------------------
3130       end module io_clust