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
 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
 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,constr_dist
1290 !      implicit none
1291 !      include 'DIMENSIONS'
1292 !      include 'sizesclu.dat'
1293 !      include 'COMMON.IOUNITS'
1294 !      include 'COMMON.TIME1'
1295 !      include 'COMMON.SBRIDGE'
1296 !      include 'COMMON.CONTROL'
1297 !      include 'COMMON.CLUSTER'
1298 !      include 'COMMON.CHAIN'
1299 !      include 'COMMON.HEADER'
1300 !      include 'COMMON.FFIELD'
1301 !      include 'COMMON.FREE'
1302 !      include 'COMMON.INTERACT'
1303       character(len=320) :: controlcard !,ucase
1304 !#ifdef MPL
1305 !      include 'COMMON.INFO'
1306 !#endif
1307       integer :: i
1308
1309       read (INP,'(a80)') titel
1310       call card_concat(controlcard,.true.)
1311
1312       call readi(controlcard,'NRES',nres_molec(1),0)
1313       call readi(controlcard,'NRES_NUCL',nres_molec(2),0)
1314       call readi(controlcard,'NRES_CAT',nres_molec(5),0)
1315       nres=0
1316       do i=1,5
1317        nres=nres_molec(i)+nres
1318       enddo
1319       print *,"TU",nres_molec(:)
1320
1321 !      call alloc_clust_arrays
1322       allocate(rcutoff(max_cut+1)) !(max_cut+1)
1323       allocate(beta_h(maxT)) !(maxT)
1324       allocate(mult(nres)) !(maxres)
1325
1326
1327       call readi(controlcard,'RESCALE',rescale_mode,2)
1328       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1329       write (iout,*) "DISTCHAINMAX",distchainmax
1330       call readi(controlcard,'PDBOUT',outpdb,0)
1331       call readi(controlcard,'MOL2OUT',outmol2,0)
1332       refstr=(index(controlcard,'REFSTR').gt.0)
1333       write (iout,*) "REFSTR",refstr
1334       pdbref=(index(controlcard,'PDBREF').gt.0)
1335       iscode=index(controlcard,'ONE_LETTER')
1336       tree=(index(controlcard,'MAKE_TREE').gt.0)
1337       min_var=(index(controlcard,'MINVAR').gt.0)
1338       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1339       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1340       call readi(controlcard,'NCUT',ncut,1)
1341       call readi(controlcard,'SYM',symetr,1)
1342       write (iout,*) 'sym', symetr
1343       call readi(controlcard,'NSTART',nstart,0)
1344       call reada(controlcard,'BOXX',boxxsize,100.0d0)
1345       call reada(controlcard,'BOXY',boxysize,100.0d0)
1346       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1347       call readi(controlcard,'NCLUST',nclust,5)
1348 !      ions=index(controlcard,"IONS").gt.0
1349       call readi(controlcard,"SCELEMODE",scelemode,0)
1350       print *,"SCELE",scelemode
1351       call readi(controlcard,'TORMODE',tor_mode,0)
1352 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1353         write(iout,*) "torsional and valence angle mode",tor_mode
1354       call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
1355       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1356       write(iout,*) "R_CUT_ELE=",r_cut_ele
1357       call readi(controlcard,'NEND',nend,0)
1358       call reada(controlcard,'ECUT',ecut,10.0d0)
1359       call reada(controlcard,'PROB',prob_limit,0.99d0)
1360       write (iout,*) "Probability limit",prob_limit
1361       lgrp=(index(controlcard,'LGRP').gt.0)
1362       caonly=(index(controlcard,'CA_ONLY').gt.0)
1363       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1364       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1365       call readi(controlcard,'IOPT',iopt,2)
1366       lside = index(controlcard,"SIDE").gt.0
1367       efree = index(controlcard,"EFREE").gt.0
1368       call readi(controlcard,'NTEMP',nT,1)
1369       write (iout,*) "nT",nT
1370 !el      call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1371       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1372       write (iout,*) "nT",nT
1373       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1374       do i=1,nT
1375         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1376       enddo
1377       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1378       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1379       lprint_int=index(controlcard,"PRINT_INT") .gt.0
1380       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
1381       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
1382 !c      if (constr_homology) tole=dmax1(tole,1.5d0)
1383       write (iout,*) "with_homology_constr ",with_dihed_constr,&
1384       " CONSTR_HOMOLOGY",constr_homology
1385       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
1386       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
1387       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
1388       write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
1389       if (min_var) iopt=1
1390       return
1391       end subroutine read_control
1392 !-----------------------------------------------------------------------------
1393       subroutine read_constr_homology
1394       use MPI_data
1395       use energy_data
1396       use geometry_data
1397       use control_data
1398       use control, only:init_int_table,homology_partition
1399       use MD_data, only:iset
1400 !      implicit none
1401 !      include 'DIMENSIONS'
1402 !#ifdef MPI
1403 !      include 'mpif.h'
1404 !#endif
1405 !      include 'COMMON.SETUP'
1406 !      include 'COMMON.CONTROL'
1407 !      include 'COMMON.HOMOLOGY'
1408 !      include 'COMMON.CHAIN'
1409 !      include 'COMMON.IOUNITS'
1410 !      include 'COMMON.MD'
1411 !      include 'COMMON.QRESTR'
1412 !      include 'COMMON.GEO'
1413 !      include 'COMMON.INTERACT'
1414 !      include 'COMMON.NAMES'
1415 !      include 'COMMON.VAR'
1416 !
1417
1418 !     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1419 !    &                 dist_cut
1420 !     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1421 !    &    sigma_odl_temp(maxres,maxres,max_template)
1422       character*2 kic2
1423       character*24 model_ki_dist, model_ki_angle
1424       character*500 controlcard
1425       integer :: ki,i,ii,j,k,l
1426       integer, dimension (:), allocatable :: ii_in_use
1427       integer :: i_tmp,idomain_tmp,&
1428       irec,ik,iistart,nres_temp
1429 !      integer :: iset
1430 !      external :: ilen
1431       logical :: liiflag,lfirst
1432       integer :: i01,i10
1433 !
1434 !     FP - Nov. 2014 Temporary specifications for new vars
1435 !
1436       real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1437                        rescore3_tmp, dist_cut
1438       real(kind=8), dimension (:,:),allocatable :: rescore
1439       real(kind=8), dimension (:,:),allocatable :: rescore2
1440       real(kind=8), dimension (:,:),allocatable :: rescore3
1441       real(kind=8) :: distal
1442       character*24 tpl_k_rescore
1443       character*256 pdbfile
1444
1445 ! -----------------------------------------------------------------
1446 ! Reading multiple PDB ref structures and calculation of retraints
1447 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1448 ! FP (Nov., 2014)
1449 ! -----------------------------------------------------------------
1450 !
1451 !
1452 ! Alternative: reading from input
1453       call card_concat(controlcard,.true.)
1454       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1455       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1456       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1457       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1458       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1459       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1460       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1461       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1462       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1463       if(.not.read2sigma.and.start_from_model) then
1464           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1465            write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1466           start_from_model=.false.
1467           iranconf=(indpdb.le.0)
1468       endif
1469       if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1470          write(iout,*) 'START_FROM_MODELS is ON'
1471 !      if(start_from_model .and. rest) then
1472 !        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1473 !         write(iout,*) 'START_FROM_MODELS is OFF'
1474 !         write(iout,*) 'remove restart keyword from input'
1475 !        endif
1476 !      endif
1477       if (start_from_model) nmodel_start=constr_homology
1478       if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1479       if (homol_nset.gt.1)then
1480          call card_concat(controlcard,.true.)
1481          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1482          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1483 !          write(iout,*) "iset homology_weight "
1484           do i=1,homol_nset
1485            write(iout,*) i,waga_homology(i)
1486           enddo
1487          endif
1488          iset=mod(kolor,homol_nset)+1
1489       else
1490        iset=1
1491        waga_homology(1)=1.0
1492       endif
1493
1494 !d      write (iout,*) "nnt",nnt," nct",nct
1495 !d      call flush(iout)
1496
1497       if (read_homol_frag) then
1498         call read_klapaucjusz
1499       else
1500
1501       lim_odl=0
1502       lim_dih=0
1503 !
1504 !      write(iout,*) 'nnt=',nnt,'nct=',nct
1505 !
1506 !      do i = nnt,nct
1507 !        do k=1,constr_homology
1508 !         idomain(k,i)=0
1509 !        enddo
1510 !      enddo
1511        idomain=0
1512
1513 !      ii=0
1514 !      do i = nnt,nct-2
1515 !        do j=i+2,nct
1516 !        ii=ii+1
1517 !        ii_in_use(ii)=0
1518 !        enddo
1519 !      enddo
1520       ii_in_use=0
1521       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1522       if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
1523
1524       do k=1,constr_homology
1525
1526         read(inp,'(a)') pdbfile
1527         pdbfiles_chomo(k)=pdbfile
1528         if(me.eq.king .or. .not. out1file) &
1529          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1530         pdbfile(:ilen(pdbfile))
1531         open(ipdbin,file=pdbfile,status='old',err=33)
1532         goto 34
1533   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
1534         pdbfile(:ilen(pdbfile))
1535         stop
1536   34    continue
1537 !        print *,'Begin reading pdb data'
1538 !
1539 ! Files containing res sim or local scores (former containing sigmas)
1540 !
1541
1542         write(kic2,'(bz,i2.2)') k
1543
1544         tpl_k_rescore="template"//kic2//".sco"
1545         write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1546         unres_pdb=.false.
1547         nres_temp=nres
1548         write(iout,*) "read2sigma",read2sigma
1549
1550         if (read2sigma) then
1551           call readpdb_template(k)
1552         else
1553           call readpdb
1554         endif
1555         write(iout,*) "after readpdb"
1556         if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1557         nres_chomo(k)=nres
1558         nres=nres_temp
1559         if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1560         if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1561         if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1562         if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1563         if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1564         if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1565         if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1566         if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1567         if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1568         if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1569         if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1570         if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1571         if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1572         if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1573 !        if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1574         if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1575         if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1576         if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1577         if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1578 !        if(.not.allocated(distance)) allocate(distance(constr_homology))
1579 !        if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1580
1581
1582 !
1583 !     Distance restraints
1584 !
1585 !          ... --> odl(k,ii)
1586 ! Copy the coordinates from reference coordinates (?)
1587         do i=1,2*nres_chomo(k)
1588           do j=1,3
1589             c(j,i)=cref(j,i,1)
1590 !           write (iout,*) "c(",j,i,") =",c(j,i)
1591           enddo
1592         enddo
1593 !
1594 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1595 !
1596 !         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1597           open (ientin,file=tpl_k_rescore,status='old')
1598           if (nnt.gt.1) rescore(k,1)=0.0d0
1599           do irec=nnt,nct ! loop for reading res sim
1600             if (read2sigma) then
1601              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1602                                      rescore3_tmp,idomain_tmp
1603              i_tmp=i_tmp+nnt-1
1604              idomain(k,i_tmp)=idomain_tmp
1605              rescore(k,i_tmp)=rescore_tmp
1606              rescore2(k,i_tmp)=rescore2_tmp
1607              rescore3(k,i_tmp)=rescore3_tmp
1608              if (.not. out1file .or. me.eq.king)&
1609              write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
1610                            i_tmp,rescore2_tmp,rescore_tmp,&
1611                                      rescore3_tmp,idomain_tmp
1612             else
1613              idomain(k,irec)=1
1614              read (ientin,*,end=1401) rescore_tmp
1615
1616 !           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1617              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1618 !           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1619             endif
1620           enddo
1621  1401   continue
1622         close (ientin)
1623         if (waga_dist.ne.0.0d0) then
1624           ii=0
1625           do i = nnt,nct-2
1626             do j=i+2,nct
1627
1628               x12=c(1,i)-c(1,j)
1629               y12=c(2,i)-c(2,j)
1630               z12=c(3,i)-c(3,j)
1631               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1632 !              write (iout,*) k,i,j,distal,dist2_cut
1633
1634             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
1635                  .and. distal.le.dist2_cut ) then
1636
1637               ii=ii+1
1638               ii_in_use(ii)=1
1639               l_homo(k,ii)=.true.
1640
1641 !             write (iout,*) "k",k
1642 !             write (iout,*) "i",i," j",j," constr_homology",
1643 !    &                       constr_homology
1644               ires_homo(ii)=i
1645               jres_homo(ii)=j
1646               odl(k,ii)=distal
1647               if (read2sigma) then
1648                 sigma_odl(k,ii)=0
1649                 do ik=i,j
1650                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1651                 enddo
1652                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1653                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
1654               sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1655               else
1656                 if (odl(k,ii).le.dist_cut) then
1657                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1658                 else
1659 #ifdef OLDSIGMA
1660                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
1661                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1662 #else
1663                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
1664                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1665 #endif
1666                 endif
1667               endif
1668               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1669             else
1670 !              ii=ii+1
1671 !              l_homo(k,ii)=.false.
1672             endif
1673             enddo
1674           enddo
1675         lim_odl=ii
1676         endif
1677 !        write (iout,*) "Distance restraints set"
1678 !        call flush(iout)
1679 !
1680 !     Theta, dihedral and SC retraints
1681 !
1682         if (waga_angle.gt.0.0d0) then
1683 !         open (ientin,file=tpl_k_sigma_dih,status='old')
1684 !         do irec=1,maxres-3 ! loop for reading sigma_dih
1685 !            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1686 !            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1687 !            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1688 !    &                            sigma_dih(k,i+nnt-1)
1689 !         enddo
1690 !1402   continue
1691 !         close (ientin)
1692           do i = nnt+3,nct
1693             if (idomain(k,i).eq.0) then
1694                sigma_dih(k,i)=0.0
1695                cycle
1696             endif
1697             dih(k,i)=phiref(i) ! right?
1698 !           read (ientin,*) sigma_dih(k,i) ! original variant
1699 !             write (iout,*) "dih(",k,i,") =",dih(k,i)
1700 !             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1701 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1702 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1703 !    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1704
1705             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
1706                           rescore(k,i-2)+rescore(k,i-3))/4.0
1707 !            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1708 !           write (iout,*) "Raw sigmas for dihedral angle restraints"
1709 !           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1710 !           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1711 !                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1712 !   Instead of res sim other local measure of b/b str reliability possible
1713             if (sigma_dih(k,i).ne.0) &
1714             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1715 !           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1716           enddo
1717           lim_dih=nct-nnt-2
1718         endif
1719 !        write (iout,*) "Dihedral angle restraints set"
1720 !        call flush(iout)
1721
1722         if (waga_theta.gt.0.0d0) then
1723 !         open (ientin,file=tpl_k_sigma_theta,status='old')
1724 !         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1725 !            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1726 !            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1727 !    &                              sigma_theta(k,i+nnt-1)
1728 !         enddo
1729 !1403   continue
1730 !         close (ientin)
1731
1732           do i = nnt+2,nct ! right? without parallel.
1733 !         do i = i=1,nres ! alternative for bounds acc to readpdb?
1734 !         do i=ithet_start,ithet_end ! with FG parallel.
1735              if (idomain(k,i).eq.0) then
1736               sigma_theta(k,i)=0.0
1737               cycle
1738              endif
1739              thetatpl(k,i)=thetaref(i)
1740 !            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1741 !            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1742 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1743 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1744 !            read (ientin,*) sigma_theta(k,i) ! 1st variant
1745              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
1746                              rescore(k,i-2))/3.0
1747 !             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1748              if (sigma_theta(k,i).ne.0) &
1749              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1750
1751 !            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1752 !                             rescore(k,i-2) !  right expression ?
1753 !            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1754           enddo
1755         endif
1756 !        write (iout,*) "Angle restraints set"
1757 !        call flush(iout)
1758
1759         if (waga_d.gt.0.0d0) then
1760 !       open (ientin,file=tpl_k_sigma_d,status='old')
1761 !         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1762 !            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1763 !            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1764 !    &                          sigma_d(k,i+nnt-1)
1765 !         enddo
1766 !1404   continue
1767
1768           do i = nnt,nct ! right? without parallel.
1769 !         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1770 !         do i=loc_start,loc_end ! with FG parallel.
1771                if (itype(i,1).eq.10) cycle
1772                if (idomain(k,i).eq.0 ) then
1773                   sigma_d(k,i)=0.0
1774                   cycle
1775                endif
1776                xxtpl(k,i)=xxref(i)
1777                yytpl(k,i)=yyref(i)
1778                zztpl(k,i)=zzref(i)
1779 !              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1780 !              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1781 !              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1782 !              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1783                sigma_d(k,i)=rescore3(k,i) !  right expression ?
1784                if (sigma_d(k,i).ne.0) &
1785                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1786
1787 !              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1788 !              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1789 !              read (ientin,*) sigma_d(k,i) ! 1st variant
1790           enddo
1791         endif
1792       enddo
1793 !      write (iout,*) "SC restraints set"
1794 !      call flush(iout)
1795 !
1796 ! remove distance restraints not used in any model from the list
1797 ! shift data in all arrays
1798 !
1799 !      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
1800       if (waga_dist.ne.0.0d0) then
1801         ii=0
1802         liiflag=.true.
1803         lfirst=.true.
1804         do i=nnt,nct-2
1805          do j=i+2,nct
1806           ii=ii+1
1807 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1808 !     &            .and. distal.le.dist2_cut ) then
1809 !          write (iout,*) "i",i," j",j," ii",ii
1810 !          call flush(iout)
1811           if (ii_in_use(ii).eq.0.and.liiflag.or. &
1812           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
1813             liiflag=.false.
1814             i10=ii
1815             if (lfirst) then
1816               lfirst=.false.
1817               iistart=ii
1818             else
1819               if(i10.eq.lim_odl) i10=i10+1
1820               do ki=0,i10-i01-1
1821                ires_homo(iistart+ki)=ires_homo(ki+i01)
1822                jres_homo(iistart+ki)=jres_homo(ki+i01)
1823                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
1824                do k=1,constr_homology
1825                 odl(k,iistart+ki)=odl(k,ki+i01)
1826                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
1827                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
1828                enddo
1829               enddo
1830               iistart=iistart+i10-i01
1831             endif
1832           endif
1833           if (ii_in_use(ii).ne.0.and..not.liiflag) then
1834              i01=ii
1835              liiflag=.true.
1836           endif
1837          enddo
1838         enddo
1839         lim_odl=iistart-1
1840       endif
1841 !      write (iout,*) "Removing distances completed"
1842 !      call flush(iout)
1843       endif ! .not. klapaucjusz
1844
1845       if (constr_homology.gt.0) call homology_partition
1846       write (iout,*) "After homology_partition"
1847       call flush(iout)
1848       if (constr_homology.gt.0) call init_int_table
1849       write (iout,*) "After init_int_table"
1850       call flush(iout)
1851 !      endif ! .not. klapaucjusz
1852 !      endif
1853 !      if (constr_homology.gt.0) call homology_partition
1854 !      write (iout,*) "After homology_partition"
1855 !      call flush(iout)
1856 !      if (constr_homology.gt.0) call init_int_table
1857 !      write (iout,*) "After init_int_table"
1858 !      call flush(iout)
1859 !      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1860 !      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1861 !
1862 ! Print restraints
1863 !
1864       if (.not.out_template_restr) return
1865 !d      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1866       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1867        write (iout,*) "Distance restraints from templates"
1868        do ii=1,lim_odl
1869        write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
1870         ii,ires_homo(ii),jres_homo(ii),&
1871         (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
1872         ki=1,constr_homology)
1873        enddo
1874        write (iout,*) "Dihedral angle restraints from templates"
1875        do i=nnt+3,nct
1876         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
1877             (rad2deg*dih(ki,i),&
1878             rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1879        enddo
1880        write (iout,*) "Virtual-bond angle restraints from templates"
1881        do i=nnt+2,nct
1882         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
1883             (rad2deg*thetatpl(ki,i),&
1884             rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1885        enddo
1886        write (iout,*) "SC restraints from templates"
1887        do i=nnt,nct
1888         write(iout,'(i7,100(4f8.2,4x))') i,&
1889         (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
1890          1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1891        enddo
1892       endif
1893       return
1894       end subroutine read_constr_homology
1895 !----------------------------------------------------------------------------
1896       subroutine molread
1897 !
1898 ! Read molecular data.
1899 !
1900       use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc,&
1901                   deg2rad,phibound,crefjlee,cref,rad2deg
1902       use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1903 !                 wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1904 !                 wturn3,wturn4,wturn6,wvdwpp,weights
1905       use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1906                  indpdb,constr_dist,raw_psipred, with_theta_constr
1907       use geometry, only: chainbuild,alloc_geo_arrays
1908       use energy, only: alloc_ener_arrays
1909       use io_base, only: rescode
1910       use control, only: setup_var,init_int_table
1911       use conform_compar, only: contact
1912 !      implicit none
1913 !      include 'DIMENSIONS'
1914 !      include 'COMMON.IOUNITS'
1915 !      include 'COMMON.GEO'
1916 !      include 'COMMON.VAR'
1917 !      include 'COMMON.INTERACT'
1918 !      include 'COMMON.LOCAL'
1919 !      include 'COMMON.NAMES'
1920 !      include 'COMMON.CHAIN'
1921 !      include 'COMMON.FFIELD'
1922 !      include 'COMMON.SBRIDGE'
1923 !      include 'COMMON.HEADER'
1924 !      include 'COMMON.CONTROL'
1925 !      include 'COMMON.CONTACTS'
1926 !      include 'COMMON.TIME1'
1927 !#ifdef MPL
1928 !      include 'COMMON.INFO'
1929 !#endif
1930       character(len=4) ::  sequence(nres) !(maxres)
1931       character(len=800) :: weightcard,controlcard
1932 !      integer rescode
1933       real(kind=8) :: x(6*nres) !(maxvar)
1934       real(kind=8) :: ssscale
1935       real(kind=8) :: phihel,phibet,sigmahel,sigmabet,sumv,&
1936                       secprob(3,maxres)
1937       integer  :: itype_pdb(nres) !(maxres)
1938 !      logical seq_comp
1939       integer :: i,j,kkk,mnum
1940 !
1941 ! Body
1942 !
1943 !el      allocate(weights(n_ene))
1944       allocate(weights(max_ene))
1945       call alloc_geo_arrays
1946       call alloc_ener_arrays
1947 !-----------------------------
1948       allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1949       allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1950       allocate(itype(nres+2,5)) !(maxres)
1951       allocate(molnum(nres+1))
1952       allocate(itel(nres+2))
1953
1954       do i=1,2*nres+2
1955         do j=1,3
1956           c(j,i)=0
1957           dc(j,i)=0
1958         enddo
1959       enddo
1960       do j=1,5
1961       do i=1,nres+2
1962         itype(i,j)=0
1963         itel(i)=0
1964       enddo
1965       enddo
1966 !--------------------------
1967 ! Read weights of the subsequent energy terms.
1968       call card_concat(weightcard,.true.)
1969       call reada(weightcard,'WSC',wsc,1.0d0)
1970       call reada(weightcard,'WLONG',wsc,wsc)
1971       call reada(weightcard,'WSCP',wscp,1.0d0)
1972       call reada(weightcard,'WELEC',welec,1.0D0)
1973       call reada(weightcard,'WVDWPP',wvdwpp,welec)
1974       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1975       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1976       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1977       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1978       call reada(weightcard,'WTURN3',wturn3,1.0D0)
1979       call reada(weightcard,'WTURN4',wturn4,1.0D0)
1980       call reada(weightcard,'WTURN6',wturn6,1.0D0)
1981       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1982       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1983       call reada(weightcard,'WBOND',wbond,1.0D0)
1984       call reada(weightcard,'WTOR',wtor,1.0D0)
1985       call reada(weightcard,'WTORD',wtor_d,1.0D0)
1986       call reada(weightcard,'WANG',wang,1.0D0)
1987       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1988       call reada(weightcard,'SCAL14',scal14,0.4D0)
1989       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1990       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1991       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1992       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
1993       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
1994       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
1995       if (index(weightcard,'SOFT').gt.0) ipot=6
1996       call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
1997       call reada(weightcard,'WELPP',welpp,0.0d0)
1998       call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
1999       call reada(weightcard,'WELPSB',welpsb,0.0D0)
2000       call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
2001       call reada(weightcard,'WELSB',welsb,0.0D0)
2002       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
2003       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
2004       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
2005       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
2006 !      print *,"KUR..",wtor_nucl
2007       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
2008       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
2009       call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
2010 ! 12/1/95 Added weight for the multi-body term WCORR
2011       call reada(weightcard,'WCORRH',wcorr,1.0D0)
2012       call reada(weightcard,'WSCBASE',wscbase,0.0D0)
2013       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
2014       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
2015       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
2016       call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
2017
2018       call reada(weightcard,"D0CM",d0cm,3.78d0)
2019       call reada(weightcard,"AKCM",akcm,15.1d0)
2020       call reada(weightcard,"AKTH",akth,11.0d0)
2021       call reada(weightcard,"AKCT",akct,12.0d0)
2022       call reada(weightcard,"V1SS",v1ss,-1.08d0)
2023       call reada(weightcard,"V2SS",v2ss,7.61d0)
2024       call reada(weightcard,"V3SS",v3ss,13.7d0)
2025       call reada(weightcard,"EBR",ebr,-5.50D0)
2026       call reada(weightcard,"ATRISS",atriss,0.301D0)
2027       call reada(weightcard,"BTRISS",btriss,0.021D0)
2028       call reada(weightcard,"CTRISS",ctriss,1.001D0)
2029       call reada(weightcard,"DTRISS",dtriss,1.001D0)
2030       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
2031       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
2032
2033       call reada(weightcard,"HT",Ht,0.0D0)
2034       if (dyn_ss) then
2035        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
2036         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
2037         akcm=akcm/wsc*ssscale
2038         akth=akth/wsc*ssscale
2039         akct=akct/wsc*ssscale
2040         v1ss=v1ss/wsc*ssscale
2041         v2ss=v2ss/wsc*ssscale
2042         v3ss=v3ss/wsc*ssscale
2043       else
2044         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
2045       endif
2046
2047       if (wcorr4.gt.0.0d0) wcorr=wcorr4
2048       weights(1)=wsc
2049       weights(2)=wscp
2050       weights(3)=welec
2051       weights(4)=wcorr
2052       weights(5)=wcorr5
2053       weights(6)=wcorr6
2054       weights(7)=wel_loc
2055       weights(8)=wturn3
2056       weights(9)=wturn4
2057       weights(10)=wturn6
2058       weights(11)=wang
2059       weights(12)=wscloc
2060       weights(13)=wtor
2061       weights(14)=wtor_d
2062       weights(15)=wstrain
2063       weights(16)=wvdwpp
2064       weights(17)=wbond
2065       weights(18)=scal14
2066 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
2067       weights(21)=wsccor
2068           weights(26)=wvdwpp_nucl
2069           weights(27)=welpp
2070           weights(28)=wvdwpsb
2071           weights(29)=welpsb
2072           weights(30)=wvdwsb
2073           weights(31)=welsb
2074           weights(32)=wbond_nucl
2075           weights(33)=wang_nucl
2076           weights(34)=wsbloc
2077           weights(35)=wtor_nucl
2078           weights(36)=wtor_d_nucl
2079           weights(37)=wcorr_nucl
2080           weights(38)=wcorr3_nucl
2081           weights(41)=wcatcat
2082           weights(42)=wcatprot
2083           weights(46)=wscbase
2084           weights(47)=wscpho
2085           weights(48)=wpeppho
2086           weights(49)=wpeppho
2087           weights(50)=wcatnucl
2088
2089
2090       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
2091         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
2092         wturn4,wturn6,wsccor
2093    10 format (/'Energy-term weights (unscaled):'// &
2094        'WSCC=   ',f10.6,' (SC-SC)'/ &
2095        'WSCP=   ',f10.6,' (SC-p)'/ &
2096        'WELEC=  ',f10.6,' (p-p electr)'/ &
2097        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
2098        'WBOND=  ',f10.6,' (stretching)'/ &
2099        'WANG=   ',f10.6,' (bending)'/ &
2100        'WSCLOC= ',f10.6,' (SC local)'/ &
2101        'WTOR=   ',f10.6,' (torsional)'/ &
2102        'WTORD=  ',f10.6,' (double torsional)'/ &
2103        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
2104        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
2105        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
2106        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
2107        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
2108        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
2109        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
2110        'WTURN6= ',f10.6,' (turns, 6th order)'/ &
2111        'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
2112
2113       if (wcorr4.gt.0.0d0) then
2114         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
2115          'between contact pairs of peptide groups'
2116         write (iout,'(2(a,f5.3/))') &
2117         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
2118         'Range of quenching the correlation terms:',2*delt_corr
2119       else if (wcorr.gt.0.0d0) then
2120         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
2121          'between contact pairs of peptide groups'
2122       endif
2123       write (iout,'(a,f8.3)') &
2124         'Scaling factor of 1,4 SC-p interactions:',scal14
2125       write (iout,'(a,f8.3)') &
2126         'General scaling factor of SC-p interactions:',scalscp
2127       r0_corr=cutoff_corr-delt_corr
2128       do i=1,20
2129         aad(i,1)=scalscp*aad(i,1)
2130         aad(i,2)=scalscp*aad(i,2)
2131         bad(i,1)=scalscp*bad(i,1)
2132         bad(i,2)=scalscp*bad(i,2)
2133       enddo
2134
2135       call flush(iout)
2136       print *,'indpdb=',indpdb,' pdbref=',pdbref
2137
2138 ! Read sequence if not taken from the pdb file.
2139       if (iscode.gt.0) then
2140         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
2141       else
2142         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
2143       endif
2144       print *,nres_molec(:),nres
2145 ! Convert sequence to numeric code
2146       do i=1,nres_molec(1)
2147         mnum=1
2148         molnum(i)=1
2149         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2150       enddo
2151       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
2152         mnum=2
2153         molnum(i)=2
2154         write (iout,*),i,sequence(i)
2155         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2156
2157       enddo
2158       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
2159         mnum=5
2160         molnum(i)=5
2161         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2162       enddo
2163       print *,nres
2164       print '(20i4)',(itype(i,molnum(i)),i=1,nres)
2165
2166       do i=1,nres-1
2167 #ifdef PROCOR
2168         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
2169             .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
2170 #else
2171         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
2172 #endif
2173           itel(i)=0
2174 #ifdef PROCOR
2175         else if (iabs(itype(i+1,1)).ne.20) then
2176 #else
2177         else if (iabs(itype(i,1)).ne.20) then
2178 #endif
2179           itel(i)=1
2180         else
2181           itel(i)=2
2182         endif
2183       enddo
2184       write (iout,*) "ITEL"
2185       do i=1,nres-1
2186         write (iout,*) i,itype(i,molnum(i)),itel(i)
2187       enddo
2188
2189       print *,'Call Read_Bridge.'
2190       call read_bridge
2191       nnt=1
2192       nct=nres
2193       print *,'NNT=',NNT,' NCT=',NCT
2194       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
2195       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
2196       if (nstart.lt.nnt) nstart=nnt
2197       if (nend.gt.nct .or. nend.eq.0) nend=nct
2198       write (iout,*) "nstart",nstart," nend",nend
2199       nres0=nres
2200 !      if (pdbref) then
2201 !        read(inp,'(a)') pdbfile
2202 !        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
2203 !        open(ipdbin,file=pdbfile,status='old',err=33)
2204 !        goto 34
2205 !  33    write (iout,'(a)') 'Error opening PDB file.'
2206 !        stop
2207 !  34    continue
2208 !        print *,'Begin reading pdb data'
2209 !        call readpdb
2210 !        print *,'Finished reading pdb data'
2211 !        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
2212 !        do i=1,nres
2213 !          itype_pdb(i)=itype(i)
2214 !        enddo
2215 !        close (ipdbin)
2216 !        write (iout,'(a,i3)') 'nsup=',nsup
2217 !        nstart_seq=nnt
2218 !        if (nsup.le.(nct-nnt+1)) then
2219 !          do i=0,nct-nnt+1-nsup
2220 !            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
2221 !              nstart_seq=nnt+i
2222 !              goto 111
2223 !            endif
2224 !          enddo
2225 !          write (iout,'(a)')
2226 !     &            'Error - sequences to be superposed do not match.'
2227 !          stop
2228 !        else
2229 !          do i=0,nsup-(nct-nnt+1)
2230 !            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
2231 !     &      then
2232 !              nstart_sup=nstart_sup+i
2233 !              nsup=nct-nnt+1
2234 !              goto 111
2235 !            endif
2236 !          enddo
2237 !          write (iout,'(a)')
2238 !     &            'Error - sequences to be superposed do not match.'
2239 !        endif
2240 !  111   continue
2241 !        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
2242 !     &                 ' nstart_seq=',nstart_seq
2243 !      endif
2244             if (with_dihed_constr) then
2245
2246       read (inp,*) ndih_constr
2247       if (ndih_constr.gt.0) then
2248         raw_psipred=.false.
2249 !C        read (inp,*) ftors
2250 !C        write (iout,*) 'FTORS',ftors
2251 !C ftors is the force constant for torsional quartic constrains
2252         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),&
2253                      i=1,ndih_constr)
2254         write (iout,*)&
2255         'There are',ndih_constr,' constraints on phi angles.'
2256         do i=1,ndih_constr
2257           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),&
2258        ftors(i)
2259         enddo
2260         do i=1,ndih_constr
2261           phi0(i)=deg2rad*phi0(i)
2262           drange(i)=deg2rad*drange(i)
2263         enddo
2264       else if (ndih_constr.lt.0) then
2265         raw_psipred=.true.
2266         call card_concat(controlcard,.true.)
2267         call reada(controlcard,"PHIHEL",phihel,50.0D0)
2268         call reada(controlcard,"PHIBET",phibet,180.0D0)
2269         call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
2270         call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
2271         call reada(controlcard,"WDIHC",wdihc,0.591d0)
2272         write (iout,*) "Weight of the dihedral restraint term",wdihc
2273         read(inp,'(9x,3f7.3)')&
2274           (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
2275         write (iout,*) "The secprob array"
2276         do i=nnt,nct
2277           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
2278         enddo
2279         ndih_constr=0
2280         do i=nnt+3,nct
2281           if (itype(i-3,molnum(i-3)).ne.ntyp1 .and. itype(i-2,molnum(i-2)).ne.ntyp1&
2282          .and. itype(i-1,molnum(i-1)).ne.ntyp1 .and. itype(i,molnum(i)).ne.ntyp1) then
2283             ndih_constr=ndih_constr+1
2284             idih_constr(ndih_constr)=i
2285             sumv=0.0d0
2286             do j=1,3
2287               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
2288               sumv=sumv+vpsipred(j,ndih_constr)
2289             enddo
2290             do j=1,3
2291               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
2292             enddo
2293             phibound(1,ndih_constr)=phihel*deg2rad
2294             phibound(2,ndih_constr)=phibet*deg2rad
2295             sdihed(1,ndih_constr)=sigmahel*deg2rad
2296             sdihed(2,ndih_constr)=sigmabet*deg2rad
2297           endif
2298         enddo
2299         write (iout,*)&
2300         'There are',ndih_constr,&
2301         ' bimodal restraints on gamma angles.'
2302         do i=1,ndih_constr
2303           write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,&
2304            restyp(itype(idih_constr(i)-2,molnum(idih_constr(i)-2)),&
2305                  molnum(idih_constr(i)-2)),idih_constr(i)-2,&
2306            restyp(itype(idih_constr(i)-1,molnum(idih_constr(i)-1)),&
2307                  molnum(idih_constr(i)-2)),idih_constr(i)-1,&
2308            phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,&
2309            phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,&
2310            (vpsipred(j,i),j=1,3)
2311         enddo
2312
2313       endif ! endif ndif_constr.gt.0
2314       endif ! with_dihed_constr
2315       if (with_theta_constr) then
2316 !C with_theta_constr is keyword allowing for occurance of theta constrains
2317       read (inp,*) ntheta_constr
2318 !C ntheta_constr is the number of theta constrains
2319       if (ntheta_constr.gt.0) then
2320 !C        read (inp,*) ftors
2321         read (inp,*) (itheta_constr(i),theta_constr0(i),&
2322        theta_drange(i),for_thet_constr(i),&
2323        i=1,ntheta_constr)
2324 !C the above code reads from 1 to ntheta_constr
2325 !C itheta_constr(i) residue i for which is theta_constr
2326 !C theta_constr0 the global minimum value
2327 !C theta_drange is range for which there is no energy penalty
2328 !C for_thet_constr is the force constant for quartic energy penalty
2329 !C E=k*x**4
2330          write (iout,*)&
2331         'There are',ntheta_constr,' constraints on phi angles.'
2332          do i=1,ntheta_constr
2333           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),&
2334          theta_drange(i),&
2335          for_thet_constr(i)
2336          enddo
2337 !C        endif
2338         do i=1,ntheta_constr
2339          theta_constr0(i)=deg2rad*theta_constr0(i)
2340           theta_drange(i)=deg2rad*theta_drange(i)
2341         enddo
2342 !C        if(me.eq.king.or..not.out1file)
2343 !C     &   write (iout,*) 'FTORS',ftors
2344 !C        do i=1,ntheta_constr
2345 !C          ii = itheta_constr(i)
2346 !C          thetabound(1,ii) = phi0(i)-drange(i)
2347 !C          thetabound(2,ii) = phi0(i)+drange(i)
2348 !C        enddo
2349       endif ! ntheta_constr.gt.0
2350       endif! with_theta_constr
2351       if (constr_homology.gt.0) then
2352 !c        write (iout,*) "About to call read_constr_homology"
2353 !c        call flush(iout)
2354         call read_constr_homology
2355 !c        write (iout,*) "Exit read_constr_homology"
2356 !c        call flush(iout)
2357         if (indpdb.gt.0 .or. pdbref) then
2358           kkk=1
2359           do i=1,2*nres
2360             do j=1,3
2361               c(j,i)=crefjlee(j,i)
2362               cref(j,i,kkk)=crefjlee(j,i)
2363             enddo
2364           enddo
2365         endif
2366 #ifdef DEBUG
2367         write (iout,*) "Array C"
2368         do i=1,nres
2369           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
2370            (c(j,i+nres),j=1,3)
2371         enddo
2372         write (iout,*) "Array Cref"
2373         do i=1,nres
2374           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),&
2375            (cref(j,i+nres),j=1,3)
2376         enddo
2377 #endif
2378 #ifdef DEBUG
2379        call int_from_cart1(.false.)
2380        call sc_loc_geom(.false.)
2381        do i=1,nres
2382          thetaref(i)=theta(i)
2383          phiref(i)=phi(i)
2384          write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
2385        enddo
2386        do i=1,nres-1
2387          do j=1,3
2388            dc(j,i)=c(j,i+1)-c(j,i)
2389            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2390          enddo
2391        enddo
2392        do i=2,nres-1
2393          do j=1,3
2394            dc(j,i+nres)=c(j,i+nres)-c(j,i)
2395            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2396          enddo
2397        enddo
2398 #endif
2399       else
2400         homol_nset=0
2401       endif
2402
2403
2404       write(iout,*)"przed ini_int_tab"
2405       call init_int_table
2406       write(iout,*)"po ini_int_tab"
2407       write(iout,*)"przed setup var"
2408       call setup_var
2409       write(iout,*)"po setup var"
2410       if (ns.gt.0.and.dyn_ss) then
2411           do i=nss+1,nhpb
2412             ihpb(i-nss)=ihpb(i)
2413             jhpb(i-nss)=jhpb(i)
2414             forcon(i-nss)=forcon(i)
2415             dhpb(i-nss)=dhpb(i)
2416           enddo
2417           nhpb=nhpb-nss
2418           nss=0
2419           link_start=1
2420           link_end=nhpb
2421 !          call hpb_partition
2422           do i=1,ns
2423             dyn_ss_mask(iss(i))=.true.
2424           enddo
2425       endif
2426
2427       write (iout,*) "molread: REFSTR",refstr
2428       if (refstr) then
2429         if (.not.pdbref) then
2430           call read_angles(inp,*38)
2431           goto 39
2432    38     write (iout,'(a)') 'Error reading reference structure.'
2433 #ifdef MPL
2434           call mp_stopall(Error_Msg)
2435 #else
2436           stop 'Error reading reference structure'
2437 #endif
2438    39     call chainbuild
2439           nstart_sup=nnt
2440           nstart_seq=nnt
2441           nsup=nct-nnt+1
2442           kkk=1
2443           do i=1,2*nres
2444             do j=1,3
2445               cref(j,i,kkk)=c(j,i)
2446             enddo
2447           enddo
2448         endif
2449         call contact(.true.,ncont_ref,icont_ref)
2450       endif
2451       return
2452       end subroutine molread
2453 !-----------------------------------------------------------------------------
2454       subroutine openunits
2455 !      implicit none
2456 !      include 'DIMENSIONS'
2457       use control_data, only: from_cx,from_bx,from_cart
2458 #ifdef MPI
2459       use MPI_data
2460       include "mpif.h"
2461       character(len=3) :: liczba
2462 !      include "COMMON.MPI"
2463 #endif
2464 !      include 'COMMON.IOUNITS'
2465 !      include 'COMMON.CONTROL'
2466       integer :: lenpre,lenpot !,ilen
2467 !      external ilen
2468       character(len=16) :: cformat,cprint
2469 !      character(len=16) ucase
2470       integer :: lenint,lenout
2471       call getenv('INPUT',prefix)
2472       call getenv('OUTPUT',prefout)
2473       call getenv('INTIN',prefintin)
2474       call getenv('COORD',cformat)
2475       call getenv('PRINTCOOR',cprint)
2476       call getenv('SCRATCHDIR',scratchdir)
2477       from_bx=.true.
2478       from_cx=.false.
2479       if (index(ucase(cformat),'CX').gt.0) then
2480         from_cx=.true.
2481         from_bx=.false.
2482       endif
2483       from_cart=.true.
2484       lenpre=ilen(prefix)
2485       lenout=ilen(prefout)
2486       lenint=ilen(prefintin)
2487 ! Get the names and open the input files
2488       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
2489 #ifdef MPI
2490       write (liczba,'(bz,i3.3)') me
2491       outname=prefout(:lenout)//'_clust.out_'//liczba
2492 #else
2493       outname=prefout(:lenout)//'_clust.out'
2494 #endif
2495       if (from_bx) then
2496         intinname=prefintin(:lenint)//'.bx'
2497       else if (from_cx) then
2498         intinname=prefintin(:lenint)//'.cx'
2499       else
2500         intinname=prefintin(:lenint)//'.int'
2501       endif
2502       rmsname=prefintin(:lenint)//'.rms'
2503       open (jplot,file=prefout(:ilen(prefout))//'.tex',&
2504              status='unknown')
2505       open (jrms,file=rmsname,status='unknown')
2506       open(iout,file=outname,status='unknown')
2507 ! Get parameter filenames and open the parameter files.
2508       call getenv('BONDPAR',bondname)
2509       open (ibond,file=bondname,status='old')
2510       call getenv('THETPAR',thetname)
2511       open (ithep,file=thetname,status='old')
2512       call getenv('ROTPAR',rotname)
2513       open (irotam,file=rotname,status='old')
2514       call getenv('TORPAR',torname)
2515       open (itorp,file=torname,status='old')
2516       call getenv('TORDPAR',tordname)
2517       open (itordp,file=tordname,status='old')
2518       call getenv('FOURIER',fouriername)
2519       open (ifourier,file=fouriername,status='old')
2520       call getenv('ELEPAR',elename)
2521       open (ielep,file=elename,status='old')
2522       call getenv('SIDEPAR',sidename)
2523       open (isidep,file=sidename,status='old')
2524       call getenv('SIDEP',sidepname)
2525       open (isidep1,file=sidepname,status="old")
2526       call getenv('SCCORPAR',sccorname)
2527       open (isccor,file=sccorname,status="old")
2528       call getenv('BONDPAR_NUCL',bondname_nucl)
2529       open (ibond_nucl,file=bondname_nucl,status='old')
2530       call getenv('THETPAR_NUCL',thetname_nucl)
2531       open (ithep_nucl,file=thetname_nucl,status='old')
2532       call getenv('ROTPAR_NUCL',rotname_nucl)
2533       open (irotam_nucl,file=rotname_nucl,status='old')
2534       call getenv('TORPAR_NUCL',torname_nucl)
2535       open (itorp_nucl,file=torname_nucl,status='old')
2536       call getenv('TORDPAR_NUCL',tordname_nucl)
2537       open (itordp_nucl,file=tordname_nucl,status='old')
2538       call getenv('SIDEPAR_NUCL',sidename_nucl)
2539       open (isidep_nucl,file=sidename_nucl,status='old')
2540       call getenv('SIDEPAR_SCBASE',sidename_scbase)
2541       open (isidep_scbase,file=sidename_scbase,status='old')
2542       call getenv('PEPPAR_PEPBASE',pepname_pepbase)
2543       open (isidep_pepbase,file=pepname_pepbase,status='old')
2544       call getenv('SCPAR_PHOSPH',pepname_scpho)
2545       open (isidep_scpho,file=pepname_scpho,status='old')
2546       call getenv('PEPPAR_PHOSPH',pepname_peppho)
2547       open (isidep_peppho,file=pepname_peppho,status='old')
2548
2549
2550       call getenv('LIPTRANPAR',liptranname)
2551       open (iliptranpar,file=liptranname,status='old')
2552       call getenv('TUBEPAR',tubename)
2553       open (itube,file=tubename,status='old')
2554       call getenv('IONPAR',ionname)
2555       open (iion,file=ionname,status='old')
2556       call getenv('IONPAR_NUCL',ionnuclname)
2557       open (iionnucl,file=ionnuclname,status='old')
2558
2559 #ifndef OLDSCP
2560 !
2561 ! 8/9/01 In the newest version SCp interaction constants are read from a file
2562 ! Use -DOLDSCP to use hard-coded constants instead.
2563 !
2564       call getenv('SCPPAR',scpname)
2565       open (iscpp,file=scpname,status='old')
2566       call getenv('SCPPAR_NUCL',scpname_nucl)
2567       open (iscpp_nucl,file=scpname_nucl,status='old')
2568
2569 #endif
2570       return
2571       end subroutine openunits
2572 !-----------------------------------------------------------------------------
2573 ! geomout.F
2574 !-----------------------------------------------------------------------------
2575       subroutine pdboutC(etot,rmsd,tytul)
2576
2577       use energy_data, only: ihpb,jhpb,itype,molnum
2578       use energy, only:boxshift,to_box
2579 !      implicit real*8 (a-h,o-z)
2580 !      include 'DIMENSIONS'
2581 !      include 'COMMON.CONTROL'
2582 !      include 'COMMON.CHAIN'
2583 !      include 'COMMON.INTERACT'
2584 !      include 'COMMON.NAMES'
2585 !      include 'COMMON.IOUNITS'
2586 !      include 'COMMON.HEADER'
2587 !      include 'COMMON.SBRIDGE'
2588 !      include 'COMMON.TEMPFAC'
2589       character(len=50) :: tytul
2590       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
2591                                                'G','H','I','J'/)
2592       integer :: ica(nres),k,iti1
2593       real(kind=8) :: etot,rmsd
2594       integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
2595       real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
2596       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
2597         ' ENERGY ',etot,' RMS ',rmsd
2598       iatom=0
2599       ichain=1
2600       ires=0
2601 !      boxxxshift(1)=int(c(1,nnt)/boxxsize)
2602 !      if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
2603 !      boxxxshift(2)=int(c(2,nnt)/boxzsize)
2604 !      if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
2605 !      boxxxshift(3)=int(c(3,nnt)/boxzsize)
2606 !      if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
2607 !      do i=1,3
2608 !        if (boxxxshift(i).gt.2) boxxxshift(i)=2
2609 !        if (boxxxshift(i).lt.-2) boxxxshift(i)=-2
2610 !
2611 !      enddo
2612       do j=1,3
2613       cbeg(j)=c(j,nnt)
2614       enddo
2615       call to_box(cbeg(1),cbeg(2),cbeg(3))
2616       boxxxx(1)=boxxsize
2617       boxxxx(2)=boxysize
2618       boxxxx(3)=boxzsize
2619
2620       do i=nnt,nct
2621         mnum=molnum(i)
2622         iti=itype(i,mnum)
2623         if (i.ne.nct) then
2624         iti1=itype(i+1,mnum)
2625         mnum1=molnum(i+1)
2626         if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle
2627         endif
2628         if (iti.eq.ntyp1_molec(mnum)) then
2629           ichain=ichain+1
2630           ires=0
2631           write (ipdb,'(a)') 'TER'
2632         else
2633         ires=ires+1
2634         if ((ires.eq.2).and.(mnum.ne.5)) then
2635          do j=1,3
2636          cbeg(j)=c(j,i-1)
2637          enddo
2638         endif
2639         iatom=iatom+1
2640         ica(i)=iatom
2641         call to_box(c(1,i),c(2,i),c(3,i))
2642
2643         DO k=1,3
2644         Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k))
2645         c(k,i)=cbeg(k)+Rdist(k)
2646         ENDDO
2647       if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then
2648         call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres))
2649         DO k=1,3
2650         Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k))
2651         c(k,i+nres)=cbeg(k)+Rdist(k)
2652         ENDDO
2653        endif
2654         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
2655            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
2656         if ((iti.ne.10).and.(mnum.ne.5)) then
2657           iatom=iatom+1
2658           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
2659             ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
2660         endif
2661         endif
2662       enddo
2663       write (ipdb,'(a)') 'TER'
2664       do i=nnt,nct-1
2665         mnum=molnum(i)
2666         mnum1=molnum(i+1)
2667         if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
2668         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2669           write (ipdb,30) ica(i),ica(i+1)
2670         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2671           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2672         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
2673           write (ipdb,30) ica(i),ica(i)+1
2674         endif
2675       enddo
2676       if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
2677         write (ipdb,30) ica(nct),ica(nct)+1
2678       endif
2679       do i=1,nss
2680         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2681       enddo
2682       write (ipdb,'(a6)') 'ENDMDL'
2683   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2684   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2685   30  FORMAT ('CONECT',8I5)
2686       return
2687       end subroutine pdboutC
2688 !-----------------------------------------------------------------------------
2689       subroutine cartout(igr,i,etot,free,rmsd,plik)
2690 !     implicit real*8 (a-h,o-z)
2691 !     include 'DIMENSIONS'
2692 !     include 'sizesclu.dat'
2693 !     include 'COMMON.IOUNITS'
2694 !     include 'COMMON.CHAIN'
2695 !     include 'COMMON.VAR'
2696 !     include 'COMMON.LOCAL'
2697 !     include 'COMMON.INTERACT'
2698 !     include 'COMMON.NAMES'
2699 !     include 'COMMON.GEO'
2700 !     include 'COMMON.CLUSTER'
2701       integer :: igr,i,j,k
2702       real(kind=8) :: etot,free,rmsd
2703       character(len=80) :: plik
2704       open (igeom,file=plik,position='append')
2705       write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
2706       write (igeom,'(i4,$)') &
2707         nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
2708       write (igeom,'(i10)') iscore(i)
2709       write (igeom,'(8f10.5)') &
2710         ((allcart(k,j,i),k=1,3),j=1,nres),&
2711         ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
2712       return
2713       end subroutine cartout
2714 !------------------------------------------------------------------------------
2715       subroutine read_klapaucjusz
2716       use energy_data
2717       use control_data, only:unres_pdb
2718       use geometry_data, only:theta,thetaref,phi,phiref,&
2719       xxref,yyref,zzref
2720       implicit none
2721 !     include 'DIMENSIONS'
2722 !#ifdef MPI
2723 !     include 'mpif.h'
2724 !#endif
2725 !     include 'COMMON.SETUP'
2726 !     include 'COMMON.CONTROL'
2727 !     include 'COMMON.HOMOLOGY'
2728 !     include 'COMMON.CHAIN'
2729 !     include 'COMMON.IOUNITS'
2730 !     include 'COMMON.MD'
2731 !     include 'COMMON.GEO'
2732 !     include 'COMMON.INTERACT'
2733 !     include 'COMMON.NAMES'
2734       character(len=256) fragfile
2735       integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2736                          ii_in_use
2737       integer, dimension(:,:), allocatable :: iresclust,inclust
2738       integer :: nclust
2739
2740       character(len=2) :: kic2
2741       character(len=24) :: model_ki_dist, model_ki_angle
2742       character(len=500) :: controlcard
2743       integer :: ki, i, j, jj,k, l, i_tmp,&
2744       idomain_tmp,&
2745       ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2746       i01,i10,nnt_chain,nct_chain
2747       real(kind=8) :: distal
2748       logical :: lprn = .true.
2749       integer :: nres_temp
2750 !      integer :: ilen
2751 !      external :: ilen
2752       logical :: liiflag,lfirst
2753
2754       real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2755       real(kind=8), dimension (:,:), allocatable  :: rescore
2756       real(kind=8), dimension (:,:), allocatable :: rescore2
2757       character(len=24) :: tpl_k_rescore
2758       character(len=256) :: pdbfile
2759
2760 !
2761 ! For new homol impl
2762 !
2763 !     include 'COMMON.VAR'
2764 !
2765 !      write (iout,*) "READ_KLAPAUCJUSZ"
2766 !      print *,"READ_KLAPAUCJUSZ"
2767 !      call flush(iout)
2768       call getenv("FRAGFILE",fragfile)
2769       write (iout,*) "Opening", fragfile
2770       call flush(iout)
2771       open(ientin,file=fragfile,status="old",err=10)
2772 !      write (iout,*) " opened"
2773 !      call flush(iout)
2774
2775       sigma_theta=0.0
2776       sigma_d=0.0
2777       sigma_dih=0.0
2778       l_homo = .false.
2779
2780       nres_temp=nres
2781       itype_temp(:)=itype(:,1)
2782       ii=0
2783       lim_odl=0
2784
2785 !      write (iout,*) "Entering loop"
2786 !      call flush(iout)
2787
2788       DO IGR = 1,NCHAIN_GROUP
2789
2790 !      write (iout,*) "igr",igr
2791       call flush(iout)
2792       read(ientin,*) constr_homology,nclust
2793       if (start_from_model) then
2794         nmodel_start=constr_homology
2795       else
2796         nmodel_start=0
2797       endif
2798
2799       ii_old=lim_odl
2800
2801       ichain=iequiv(1,igr)
2802       nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2803       nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2804 !      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2805
2806 ! Read pdb files
2807       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2808       if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2809       do k=1,constr_homology
2810         read(ientin,'(a)') pdbfile
2811         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2812         pdbfile(:ilen(pdbfile))
2813         pdbfiles_chomo(k)=pdbfile
2814         open(ipdbin,file=pdbfile,status='old',err=33)
2815         goto 34
2816   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
2817         pdbfile(:ilen(pdbfile))
2818         stop
2819   34    continue
2820         unres_pdb=.false.
2821 !        nres_temp=nres
2822         call readpdb_template(k)
2823         nres_chomo(k)=nres
2824 !        nres=nres_temp
2825         do i=1,nres
2826           rescore(k,i)=0.2d0
2827           rescore2(k,i)=1.0d0
2828         enddo
2829       enddo
2830 ! Read clusters
2831       do i=1,nclust
2832         read(ientin,*) ninclust(i),nresclust(i)
2833         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2834         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2835       enddo
2836 !
2837 ! Loop over clusters
2838 !
2839       do l=1,nclust
2840         do ll = 1,ninclust(l)
2841
2842         k = inclust(ll,l)
2843 !        write (iout,*) "l",l," ll",ll," k",k
2844         do i=1,nres
2845           idomain(k,i)=0
2846         enddo
2847         do i=1,nresclust(l)
2848           if (nnt.gt.1)  then
2849             idomain(k,iresclust(i,l)+1) = 1
2850           else
2851             idomain(k,iresclust(i,l)) = 1
2852           endif
2853         enddo
2854 !
2855 !     Distance restraints
2856 !
2857 !          ... --> odl(k,ii)
2858 ! Copy the coordinates from reference coordinates (?)
2859 !        nres_temp=nres
2860         nres=nres_chomo(k)
2861         do i=1,2*nres
2862           do j=1,3
2863             c(j,i)=chomo(j,i,k)
2864 !           write (iout,*) "c(",j,i,") =",c(j,i)
2865           enddo
2866         enddo
2867         call int_from_cart(.true.,.false.)
2868         call sc_loc_geom(.false.)
2869         do i=1,nres
2870           thetaref(i)=theta(i)
2871           phiref(i)=phi(i)
2872         enddo
2873 !        nres=nres_temp
2874         if (waga_dist.ne.0.0d0) then
2875           ii=ii_old
2876 !          do i = nnt,nct-2
2877           do i = nnt_chain,nct_chain-2
2878 !            do j=i+2,nct
2879             do j=i+2,nct_chain
2880
2881               x12=c(1,i)-c(1,j)
2882               y12=c(2,i)-c(2,j)
2883               z12=c(3,i)-c(3,j)
2884               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2885 !              write (iout,*) k,i,j,distal,dist2_cut
2886
2887             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2888                  .and. distal.le.dist2_cut ) then
2889
2890               ii=ii+1
2891               ii_in_use(ii)=1
2892               l_homo(k,ii)=.true.
2893
2894 !             write (iout,*) "k",k
2895 !             write (iout,*) "i",i," j",j," constr_homology",
2896 !     &                       constr_homology
2897               ires_homo(ii)=i+chain_border1(1,igr)-1
2898               jres_homo(ii)=j+chain_border1(1,igr)-1
2899               odl(k,ii)=distal
2900               if (read2sigma) then
2901                 sigma_odl(k,ii)=0
2902                 do ik=i,j
2903                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2904                 enddo
2905                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2906                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2907              sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2908               else
2909                 if (odl(k,ii).le.dist_cut) then
2910                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2911                 else
2912 #ifdef OLDSIGMA
2913                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2914                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2915 #else
2916                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2917                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2918 #endif
2919                 endif
2920               endif
2921               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2922             else
2923               ii=ii+1
2924 !              l_homo(k,ii)=.false.
2925             endif
2926             enddo
2927           enddo
2928         lim_odl=ii
2929         endif
2930 !
2931 !     Theta, dihedral and SC retraints
2932 !
2933         if (waga_angle.gt.0.0d0) then
2934           do i = nnt_chain+3,nct_chain
2935             iii=i+chain_border1(1,igr)-1
2936             if (idomain(k,i).eq.0) then
2937 !               sigma_dih(k,i)=0.0
2938                cycle
2939             endif
2940             dih(k,iii)=phiref(i)
2941             sigma_dih(k,iii)= &
2942                (rescore(k,i)+rescore(k,i-1)+ &
2943                            rescore(k,i-2)+rescore(k,i-3))/4.0
2944 !            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2945 !     &       " sigma_dihed",sigma_dih(k,i)
2946             if (sigma_dih(k,iii).ne.0) &
2947              sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2948           enddo
2949 !          lim_dih=nct-nnt-2
2950         endif
2951
2952         if (waga_theta.gt.0.0d0) then
2953           do i = nnt_chain+2,nct_chain
2954              iii=i+chain_border1(1,igr)-1
2955              if (idomain(k,i).eq.0) then
2956 !              sigma_theta(k,i)=0.0
2957               cycle
2958              endif
2959              thetatpl(k,iii)=thetaref(i)
2960              sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
2961                               rescore(k,i-2))/3.0
2962              if (sigma_theta(k,iii).ne.0) &
2963              sigma_theta(k,iii)=1.0d0/ &
2964              (sigma_theta(k,iii)*sigma_theta(k,iii))
2965           enddo
2966         endif
2967
2968         if (waga_d.gt.0.0d0) then
2969           do i = nnt_chain,nct_chain
2970              iii=i+chain_border1(1,igr)-1
2971                if (itype(i,1).eq.10) cycle
2972                if (idomain(k,i).eq.0 ) then
2973 !                  sigma_d(k,i)=0.0
2974                   cycle
2975                endif
2976                xxtpl(k,iii)=xxref(i)
2977                yytpl(k,iii)=yyref(i)
2978                zztpl(k,iii)=zzref(i)
2979                sigma_d(k,iii)=rescore(k,i)
2980                if (sigma_d(k,iii).ne.0) &
2981                 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
2982 !               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
2983           enddo
2984         endif
2985       enddo ! l
2986       enddo ! ll
2987 !
2988 ! remove distance restraints not used in any model from the list
2989 ! shift data in all arrays
2990 !
2991 !      write (iout,*) "ii_old",ii_old
2992       if (waga_dist.ne.0.0d0) then
2993 #ifdef DEBUG
2994        write (iout,*) "Distance restraints from templates"
2995        do iii=1,lim_odl
2996        write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
2997         iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
2998         (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
2999         ki=1,constr_homology)
3000        enddo
3001 #endif
3002         ii=ii_old
3003         liiflag=.true.
3004         lfirst=.true.
3005         do i=nnt_chain,nct_chain-2
3006          do j=i+2,nct_chain
3007           ii=ii+1
3008 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3009 !     &            .and. distal.le.dist2_cut ) then
3010 !          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
3011 !          call flush(iout)
3012           if (ii_in_use(ii).eq.0.and.liiflag.or. &
3013           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3014             liiflag=.false.
3015             i10=ii
3016             if (lfirst) then
3017               lfirst=.false.
3018               iistart=ii
3019             else
3020               if(i10.eq.lim_odl) i10=i10+1
3021               do ki=0,i10-i01-1
3022                ires_homo(iistart+ki)=ires_homo(ki+i01)
3023                jres_homo(iistart+ki)=jres_homo(ki+i01)
3024                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3025                do k=1,constr_homology
3026                 odl(k,iistart+ki)=odl(k,ki+i01)
3027                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3028                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3029                enddo
3030               enddo
3031               iistart=iistart+i10-i01
3032             endif
3033           endif
3034           if (ii_in_use(ii).ne.0.and..not.liiflag) then
3035              i01=ii
3036              liiflag=.true.
3037           endif
3038          enddo
3039         enddo
3040         lim_odl=iistart-1
3041       endif
3042
3043       lll=lim_odl-ii_old
3044
3045       do i=2,nequiv(igr)
3046
3047         ichain=iequiv(i,igr)
3048
3049         do j=nnt_chain,nct_chain
3050           jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
3051           do k=1,constr_homology
3052             dih(k,jj)=dih(k,j)
3053             sigma_dih(k,jj)=sigma_dih(k,j)
3054             thetatpl(k,jj)=thetatpl(k,j)
3055             sigma_theta(k,jj)=sigma_theta(k,j)
3056             xxtpl(k,jj)=xxtpl(k,j)
3057             yytpl(k,jj)=yytpl(k,j)
3058             zztpl(k,jj)=zztpl(k,j)
3059             sigma_d(k,jj)=sigma_d(k,j)
3060           enddo
3061         enddo
3062
3063         jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
3064 !        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
3065         do j=ii_old+1,lim_odl
3066           ires_homo(j+lll)=ires_homo(j)+jj
3067           jres_homo(j+lll)=jres_homo(j)+jj
3068           do k=1,constr_homology
3069             odl(k,j+lll)=odl(k,j)
3070             sigma_odl(k,j+lll)=sigma_odl(k,j)
3071             l_homo(k,j+lll)=l_homo(k,j)
3072           enddo
3073         enddo
3074
3075         ii_old=ii_old+lll
3076         lim_odl=lim_odl+lll
3077
3078       enddo
3079
3080       ENDDO ! IGR
3081
3082       if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
3083       nres=nres_temp
3084       itype(:,1)=itype_temp(:)
3085
3086       return
3087    10 stop "Error in fragment file"
3088       end subroutine read_klapaucjusz
3089
3090 !-----------------------------------------------------------------------------
3091       end module io_clust