source/cluster/io_clust.f90

   1       module io_clust
   2 !-----------------------------------------------------------------------------
   3       use clust_data
   4       use io_units
   5 !      use names
   6       use io_base !, only: ilen
   7       use geometry_data, only: nres,c,nres_molec,boxxsize,boxysize,boxzsize
   8       use energy_data, only: nnt,nct,nss
   9       use control_data, only: lside
  10       implicit none
  11 !-----------------------------------------------------------------------------
  12 !
  13 !
  14 !-----------------------------------------------------------------------------
  15       contains
  16 !-----------------------------------------------------------------------------
  17 ! wrtclust.f
  18 !-----------------------------------------------------------------------------
  19       SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
  20
  21       use hc_, only: ioffset
  22       use control_data, only: lprint_cart,lprint_int,titel
  23       use geometry, only: int_from_cart1,nres
  24 !      implicit real*8 (a-h,o-z)
  25 !      include 'DIMENSIONS'
  26 !      include 'sizesclu.dat'
  27       integer,parameter :: num_in_line=5
  28       LOGICAL :: PRINTANG(max_cut)
  29       integer :: PRINTPDB(max_cut),printmol2(max_cut)
  30 !      include 'COMMON.CONTROL'
  31 !      include 'COMMON.HEADER'
  32 !      include 'COMMON.CHAIN'
  33 !      include 'COMMON.VAR'
  34 !      include 'COMMON.CLUSTER'
  35 !      include 'COMMON.IOUNITS'
  36 !      include 'COMMON.GEO'
  37 !      include 'COMMON.FREE'
  38 !      include 'COMMON.TEMPFAC'
  39       real(kind=8) :: rmsave(maxgr)
  40       CHARACTER(len=64) :: prefixp,NUMM,MUMM,EXTEN,extmol
  41       character(len=80) :: cfname
  42       character(len=8) :: ctemper
  43       DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,&
  44            MUMM /'000'/
  45 !      external ilen
  46       integer :: ncon,icut,ib
  47       integer :: i,ii,ii1,ii2,igr,ind1,ind2,ico,icon,&
  48                  irecord,nrecord,j,k,jj,ind,ncon_lim,ncon_out
  49       real(kind=8) :: temper,curr_dist,emin,qpart,boltz,&
  50                  ave_dim,amax_dim,emin1
  51
  52
  53       allocate(tempfac(2,nres))
  54
  55       do i=1,64
  56         cfname(i:i)=" "
  57       enddo
  58 !      print *,"calling WRTCLUST",ncon
  59 !      write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
  60       rewind 80
  61       call flush(iout)
  62       temper=1.0d0/(beta_h(ib)*1.987d-3)
  63       if (temper.lt.100.0d0) then
  64         write(ctemper,'(f3.0)') temper
  65         ctemper(3:3)=" "
  66       else if (temper.lt.1000.0) then
  67         write (ctemper,'(f4.0)') temper
  68         ctemper(4:4)=" "
  69       else
  70         write (ctemper,'(f5.0)') temper
  71         ctemper(5:5)=" "
  72       endif
  73
  74       do i=1,ncon*(ncon-1)/2
  75         read (80) diss(i)
  76       enddo
  77       close(80,status='delete')
  78 !
  79 !  PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
  80 !
  81       ii1= index(intinname,'/')
  82       ii2=ii1
  83       ii1=ii1+1
  84       do while (ii2.gt.0)
  85         ii1=ii1+ii2
  86         ii2=index(intinname(ii1:),'/')
  87       enddo
  88       ii = ii1+index(intinname(ii1:),'.')-1
  89       if (ii.eq.0) then
  90         ii=ilen(intinname)
  91       else
  92         ii=ii-1
  93       endif
  94       prefixp=intinname(ii1:ii)
  95 !d    print *,icut,printang(icut),printpdb(icut),printmol2(icut)
  96 !d    print *,'ecut=',ecut
  97       WRITE (iout,100) NGR
  98       DO 19 IGR=1,NGR
  99       WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
 100       NRECORD=LICZ(IGR)/num_in_line
 101       IND1=1
 102       DO 63 IRECORD=1,NRECORD
 103       IND2=IND1+num_in_line-1
 104       WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
 105           totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
 106       IND1=IND2+1
 107    63 CONTINUE
 108       WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
 109          totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
 110       IND1=1
 111       ICON=list_conf(NCONF(IGR,1))
 112 !      WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
 113 ! 12/8/93 Estimation of "diameters" of the subsequent families.
 114       ave_dim=0.0
 115       amax_dim=0.0
 116 !      write (iout,*) "ecut",ecut
 117       do i=2,licz(igr)
 118         ii=nconf(igr,i)
 119         if (totfree(ii)-emin .gt. ecut) goto 10
 120         do j=1,i-1
 121           jj=nconf(igr,j)
 122           if (jj.eq.1) exit
 123           if (ii.lt.jj) then
 124             ind=ioffset(ncon,ii,jj)
 125           else
 126             ind=ioffset(ncon,jj,ii)
 127           endif
 128 !          write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
 129 !     &     " ind",ind
 130           call flush(iout)
 131           curr_dist=dabs(diss(ind)+0.0d0)
 132 !          write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
 133 !     &      list_conf(jj),curr_dist
 134           if (curr_dist .gt. amax_dim) amax_dim=curr_dist
 135           ave_dim=ave_dim+curr_dist**2
 136         enddo
 137       enddo
 138    10 if (licz(igr) .gt. 1) &
 139        ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
 140       write (iout,'(/A,F8.1,A,F8.1)') &
 141        'Max. distance in the family:',amax_dim,&
 142        '; average distance in the family:',ave_dim
 143       rmsave(igr)=0.0d0
 144       qpart=0.0d0
 145       do i=1,licz(igr)
 146         icon=nconf(igr,i)
 147         boltz=dexp(-totfree(icon))
 148         rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
 149         qpart=qpart+boltz
 150       enddo
 151       rmsave(igr)=rmsave(igr)/qpart
 152       write (iout,'(a,f5.2,a)') "Average RMSD",rmsave(igr)," A"
 153    19 CONTINUE
 154       WRITE (iout,400)
 155       WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
 156 !      print *,icut,printang(icut)
 157       IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
 158         emin=totfree_gr(1)
 159 !        print *,'emin',emin,' ngr',ngr
 160         if (lprint_cart) then
 161           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 162             //"K"//".x"
 163         else
 164           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 165             //"K"//".int"
 166         endif
 167         do igr=1,ngr
 168           icon=nconf(igr,1)
 169           if (totfree_gr(igr)-emin.le.ecut) then
 170             if (lprint_cart) then
 171               call cartout(igr,icon,totfree(icon)/beta_h(ib),&
 172                 totfree_gr(igr)/beta_h(ib),&
 173                 rmstb(icon),cfname)
 174             else
 175 !              print '(a)','calling briefout'
 176               do i=1,2*nres
 177                 do j=1,3
 178                   c(j,i)=allcart(j,i,icon)
 179                 enddo
 180               enddo
 181               call int_from_cart1(.false.)
 182 !el              call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),&
 183 !el                totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),&
 184 !el                jhpb_all(1,icon),cfname)
 185               call briefout(igr,totfree(icon)/beta_h(ib),&
 186                 totfree_gr(igr))
 187 !              print '(a)','exit briefout'
 188             endif
 189           endif
 190         enddo
 191         close(igeom)
 192       ENDIF
 193       IF (PRINTPDB(ICUT).gt.0) THEN
 194 ! Write out a number of conformations from each family in PDB format and
 195 ! create InsightII command file for their displaying in different colors
 196         cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 197           //"K_"//'ave'//exten
 198         write (iout,*) "cfname",cfname
 199         OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 200         write (ipdb,'(a,f8.2)') &
 201           "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
 202         close (ipdb)
 203         I=1
 204         ICON=NCONF(1,1)
 205         EMIN=totfree_gr(I)
 206         emin1=totfree(icon)
 207         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 208 !          write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
 209           write (NUMM,'(bz,i4.4)') i
 210           ncon_lim=min0(licz(i),printpdb(icut))
 211           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 212             //"K_"//numm(:ilen(numm))//exten
 213           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 214           write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5," AVE RMSD",0pf5.2)')&
 215            i,totfree_gr(i)/beta_h(ib),rmsave(i) !'
 216 ! Write conformations of the family i to PDB files
 217           ncon_out=1
 218           do while (ncon_out.lt.printpdb(icut) .and. &
 219            ncon_out.lt.licz(i).and. &
 220            totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 221             ncon_out=ncon_out+1
 222 !            write (iout,*) i,ncon_out,nconf(i,ncon_out),
 223 !     &        totfree(nconf(i,ncon_out)),emin1,ecut
 224           enddo
 225           write (iout,*) "ncon_out",ncon_out
 226           call flush(iout)
 227           do j=1,nres
 228             tempfac(1,j)=5.0d0
 229             tempfac(2,j)=5.0d0
 230           enddo
 231           do j=1,ncon_out
 232             icon=nconf(i,j)
 233             do ii=1,2*nres
 234               do k=1,3
 235                 c(k,ii)=allcart(k,ii,icon)
 236               enddo
 237             enddo
 238             call pdboutC(totfree(icon)/beta_h(ib),rmstb(icon),titel)
 239             write (ipdb,'("TER")')
 240           enddo
 241           close(ipdb)
 242 ! Average structures and structures closest to average
 243           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 244           //"K_"//'ave'//exten
 245           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',&
 246            position="APPEND")
 247           call ave_coord(i)
 248           write (ipdb,'(a,i5)') "REMARK CLUSTER",i
 249           call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 250           write (ipdb,'("TER")')
 251           call closest_coord(i)
 252           call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 253           write (ipdb,'("TER")')
 254           close (ipdb)
 255           I=I+1
 256           ICON=NCONF(I,1)
 257           emin1=totfree(icon)
 258         ENDDO
 259       ENDIF
 260       IF (printmol2(icut).gt.0) THEN
 261 ! Write out a number of conformations from each family in PDB format and
 262 ! create InsightII command file for their displaying in different colors
 263         I=1
 264         ICON=NCONF(1,1)
 265         EMIN=ENERGY(ICON)
 266         emin1=emin
 267         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 268           write (NUMM,'(bz,i4.4)') i
 269           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
 270             //"K_"//numm(:ilen(numm))//extmol
 271           OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 272           ncon_out=1
 273           do while (ncon_out.lt.printmol2(icut) .and. &
 274            ncon_out.lt.licz(i).and. &
 275            totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 276             ncon_out=ncon_out+1
 277           enddo
 278           do j=1,ncon_out
 279             icon=nconf(i,j)
 280             do ii=1,2*nres
 281               do k=1,3
 282                 c(k,ii)=allcart(k,ii,icon)
 283               enddo
 284             enddo
 285             CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
 286           enddo
 287           CLOSE(imol2)
 288           I=I+1
 289           ICON=NCONF(I,1)
 290           emin1=totfree(icon)
 291         ENDDO
 292       ENDIF
 293   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
 294   200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,&
 295        ' CONTAINS ',I4,' CONFORMATION(S): ')
 296 ! 300 FORMAT ( 8(I4,F6.1))
 297   300 FORMAT (5(I4,1pe12.3))
 298   400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
 299   500 FORMAT (8(2I4,2X))
 300   600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
 301       RETURN
 302       END SUBROUTINE WRTCLUST
 303 !------------------------------------------------------------------------------
 304       subroutine ave_coord(igr)
 305
 306       use control_data, only:lside
 307       use regularize_, only:fitsq,matvec
 308 !      implicit none
 309 !      include 'DIMENSIONS'
 310 !      include 'sizesclu.dat'
 311 !      include 'COMMON.CONTROL'
 312 !      include 'COMMON.CLUSTER'
 313 !      include 'COMMON.CHAIN'
 314 !      include 'COMMON.INTERACT'
 315 !      include 'COMMON.VAR'
 316 !      include 'COMMON.TEMPFAC'
 317 !      include 'COMMON.IOUNITS'
 318       logical :: non_conv
 319       real(kind=8) :: przes(3),obrot(3,3)
 320       real(kind=8) :: xx(3,2*nres),yy(3,2*nres),csq(3,2*nres) !(3,maxres2)
 321       real(kind=8) :: eref
 322       integer :: i,ii,j,k,icon,jcon,igr
 323       real(kind=8) :: rms,boltz,qpart,cwork(3,2*nres),cref1(3,2*nres) !(3,maxres2)
 324 !      write (iout,*) "AVE_COORD: igr",igr
 325       jcon=nconf(igr,1)
 326       eref=totfree(jcon)
 327       boltz = dexp(-totfree(jcon)+eref)
 328       qpart=boltz
 329       do i=1,2*nres
 330         do j=1,3
 331           c(j,i)=allcart(j,i,jcon)*boltz
 332           cref1(j,i)=allcart(j,i,jcon)
 333           csq(j,i)=allcart(j,i,jcon)**2*boltz
 334         enddo
 335       enddo
 336       DO K=2,LICZ(IGR)
 337       jcon=nconf(igr,k)
 338       if (lside) then
 339         ii=0
 340         do i=nnt,nct
 341           ii=ii+1
 342           do j=1,3
 343             xx(j,ii)=allcart(j,i,jcon)
 344             yy(j,ii)=cref1(j,i)
 345           enddo
 346         enddo
 347         do i=nnt,nct
 348 !          if (itype(i).ne.10) then
 349             ii=ii+1
 350             do j=1,3
 351               xx(j,ii)=allcart(j,i+nres,jcon)
 352               yy(j,ii)=cref1(j,i+nres)
 353             enddo
 354 !          endif
 355         enddo
 356         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 357       else
 358         do i=nnt,nct
 359           do j=1,3
 360             cwork(j,i)=allcart(j,i,jcon)
 361           enddo
 362         enddo
 363         call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot &
 364              ,non_conv)
 365       endif
 366 !      write (iout,*) "rms",rms
 367 !      do i=1,3
 368 !        write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
 369 !      enddo
 370       if (rms.lt.0.0) then
 371         print *,'error, rms^2 = ',rms,icon,jcon
 372         stop
 373       endif
 374       if (non_conv) print *,non_conv,icon,jcon
 375       boltz=dexp(-totfree(jcon)+eref)
 376       qpart = qpart + boltz
 377       do i=1,2*nres
 378         do j=1,3
 379           xx(j,i)=allcart(j,i,jcon)
 380         enddo
 381       enddo
 382       call matvec(cwork,obrot,xx,2*nres)
 383       do i=1,2*nres
 384 !        write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
 385 !     &    (allcart(j,i,jcon),j=1,3)
 386         do j=1,3
 387           cwork(j,i)=cwork(j,i)+przes(j)
 388           c(j,i)=c(j,i)+cwork(j,i)*boltz
 389           csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
 390         enddo
 391       enddo
 392       ENDDO ! K
 393       do i=1,2*nres
 394         do j=1,3
 395           c(j,i)=c(j,i)/qpart
 396           csq(j,i)=csq(j,i)/qpart-c(j,i)**2
 397         enddo
 398 !        write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
 399       enddo
 400       do i=nnt,nct
 401         tempfac(1,i)=0.0d0
 402         tempfac(2,i)=0.0d0
 403         do j=1,3
 404           tempfac(1,i)=tempfac(1,i)+csq(j,i)
 405           tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
 406         enddo
 407         tempfac(1,i)=dsqrt(tempfac(1,i))
 408         tempfac(2,i)=dsqrt(tempfac(2,i))
 409       enddo
 410       return
 411       end subroutine ave_coord
 412 !------------------------------------------------------------------------------
 413       subroutine closest_coord(igr)
 414
 415       use regularize_, only:fitsq
 416 !      implicit none
 417 !      include 'DIMENSIONS'
 418 !      include 'sizesclu.dat'
 419 !      include 'COMMON.IOUNITS'
 420 !      include 'COMMON.CONTROL'
 421 !      include 'COMMON.CLUSTER'
 422 !      include 'COMMON.CHAIN'
 423 !      include 'COMMON.INTERACT'
 424 !      include 'COMMON.VAR'
 425       logical :: non_conv
 426       real(kind=8) :: przes(3),obrot(3,3)
 427       real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
 428       integer :: i,ii,j,k,icon,jcon,jconmin,igr
 429       real(kind=8) :: rms,rmsmin,cwork(3,2*nres)
 430       rmsmin=1.0d10
 431       jconmin=nconf(igr,1)
 432       DO K=1,LICZ(IGR)
 433       jcon=nconf(igr,k)
 434       if (lside) then
 435         ii=0
 436         do i=nnt,nct
 437           ii=ii+1
 438           do j=1,3
 439             xx(j,ii)=allcart(j,i,jcon)
 440             yy(j,ii)=c(j,i)
 441           enddo
 442         enddo
 443         do i=nnt,nct
 444 !          if (itype(i).ne.10) then
 445             ii=ii+1
 446             do j=1,3
 447               xx(j,ii)=allcart(j,i+nres,jcon)
 448               yy(j,ii)=c(j,i+nres)
 449             enddo
 450 !          endif
 451         enddo
 452         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 453       else
 454         do i=nnt,nct
 455           do j=1,3
 456             cwork(j,i)=allcart(j,i,jcon)
 457           enddo
 458         enddo
 459         call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot,&
 460              non_conv)
 461       endif
 462       if (rms.lt.0.0) then
 463         print *,'error, rms^2 = ',rms,icon,jcon
 464         stop
 465       endif
 466 !      write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
 467       if (non_conv) print *,non_conv,icon,jcon
 468       if (rms.lt.rmsmin) then
 469         rmsmin=rms
 470         jconmin=jcon
 471       endif
 472       ENDDO ! K
 473 !      write (iout,*) "rmsmin",rmsmin," rms",rms
 474       call flush(iout)
 475       do i=1,2*nres
 476         do j=1,3
 477           c(j,i)=allcart(j,i,jconmin)
 478         enddo
 479       enddo
 480       return
 481       end subroutine closest_coord
 482 !-----------------------------------------------------------------------------
 483 ! read_coords.F
 484 !-----------------------------------------------------------------------------
 485       subroutine read_coords(ncon,*)
 486
 487       use energy_data, only: ihpb,jhpb,max_ene
 488       use control_data, only: from_bx,from_cx
 489       use control, only: tcpu
 490 !      implicit none
 491 !      include "DIMENSIONS"
 492 !      include "sizesclu.dat"
 493 #ifdef MPI
 494       use MPI_data
 495       include "mpif.h"
 496       integer :: IERROR,ERRCODE !,STATUS(MPI_STATUS_SIZE)
 497 !      include "COMMON.MPI"
 498 #else
 499       use MPI_data, only: nprocs
 500 #endif
 501 !      include "COMMON.CONTROL"
 502 !      include "COMMON.CHAIN"
 503 !      include "COMMON.INTERACT"
 504 !      include "COMMON.IOUNITS"
 505 !      include "COMMON.VAR"
 506 !      include "COMMON.SBRIDGE"
 507 !      include "COMMON.GEO"
 508 !      include "COMMON.CLUSTER"
 509       character(len=3) :: liczba
 510       integer :: ncon
 511       integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,if,ib,&
 512         nn,nn1,inan
 513       integer :: ixdrf,iret,itmp
 514       real(kind=4) :: prec,reini,refree,rmsdev
 515       integer :: nrec,nlines,iscor,lenrec,lenrec_in
 516       real(kind=8) :: energ,t_acq !,tcpu
 517 !el      integer ilen,iroof
 518 !el      external ilen,iroof
 519       real(kind=8) :: rjunk
 520       integer :: ntot_all(0:nprocs-1) !(0:maxprocs-1)
 521       logical :: lerr
 522       real(kind=8) :: energia(0:max_ene),etot
 523       real(kind=4) :: csingle(3,2*nres+2)
 524       integer :: Previous,Next
 525       character(len=256) :: bprotfiles
 526 !      print *,"Processor",me," calls read_protein_data"
 527 #ifdef MPI
 528       if (me.eq.master) then
 529         Previous=MPI_PROC_NULL
 530       else
 531         Previous=me-1
 532       endif
 533       if (me.eq.nprocs-1) then
 534         Next=MPI_PROC_NULL
 535       else
 536         Next=me+1
 537       endif
 538 ! Set the scratchfile names
 539       write (liczba,'(bz,i3.3)') me
 540
 541       allocate(STATUS(MPI_STATUS_SIZE))
 542 #endif
 543 ! 1/27/05 AL Change stored coordinates to single precision and don't store
 544 !         energy components in the binary databases.
 545       lenrec=12*(nres+nct-nnt+1)+4*(2*nss+2)+16
 546       lenrec_in=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
 547 #ifdef DEBUG
 548       write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss", nss
 549       write (iout,*) "lenrec_in",lenrec_in
 550 #endif
 551       bprotfiles=scratchdir(:ilen(scratchdir))// &
 552              "/"//prefix(:ilen(prefix))//liczba//".xbin"
 553 ! EL
 554 ! allocate cluster arrays
 555       allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
 556       allocate(entfac(maxconf)) !(maxconf)
 557       allocate(rmstb(maxconf)) !(maxconf)
 558       allocate(allcart(3,2*nres,maxstr_proc)) !(3,maxres2,maxstr_proc)
 559       allocate(nss_all(maxstr_proc)) !(maxstr_proc)
 560       allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc)
 561       allocate(iscore(maxconf)) !(maxconf)
 562
 563
 564 #ifdef CHUJ
 565       ICON=1
 566   123 continue
 567       if (from_cart .and. .not. from_bx .and. .not. from_cx) then
 568         if (efree) then
 569         read (intin,*,end=13,err=11) energy(icon),totfree(icon),&
 570           rmstb(icon),&
 571           nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
 572           i=1,nss_all(icon)),iscore(icon)
 573         else
 574         read (intin,*,end=13,err=11) energy(icon),rmstb(icon),&
 575           nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
 576           i=1,nss_all(icon)),iscore(icon)
 577         endif
 578         read (intin,'(8f10.5)',end=13,err=10) &
 579           ((allcart(j,i,icon),j=1,3),i=1,nres),&
 580           ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
 581         print *,icon,energy(icon),nss_all(icon),rmstb(icon)
 582       else
 583         read(intin,'(a80)',end=13,err=12) lineh
 584         read(lineh(:5),*,err=8) ic
 585         if (efree) then
 586         read(lineh(6:),*,err=8) energy(icon)
 587         else
 588         read(lineh(6:),*,err=8) energy(icon)
 589         endif
 590         goto 9
 591     8   ic=1
 592         print *,'error, assuming e=1d10',lineh
 593         energy(icon)=1d10
 594         nss=0
 595     9   continue
 596 !old        read(lineh(18:),*,end=13,err=11) nss_all(icon)
 597         ii = index(lineh(15:)," ")+15
 598         read(lineh(ii:),*,end=13,err=11) nss_all(icon)
 599         IF (NSS_all(icon).LT.9) THEN
 600           read (lineh(20:),*,end=102) &
 601           (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),&
 602           iscore(icon)
 603         ELSE
 604           read (lineh(20:),*,end=102) &
 605                  (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
 606           read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),&
 607             I=9,NSS_all(icon)),iscore(icon)
 608         ENDIF
 609
 610   102   continue
 611
 612         PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
 613         call read_angles(intin,*13)
 614         do i=1,nres
 615           phiall(i,icon)=phi(i)
 616           thetall(i,icon)=theta(i)
 617           alphall(i,icon)=alph(i)
 618           omall(i,icon)=omeg(i)
 619         enddo
 620       endif
 621       ICON=ICON+1
 622       GOTO 123
 623 !
 624 ! CALCULATE DISTANCES
 625 !
 626    10 print *,'something wrong with angles'
 627       goto 13
 628    11 print *,'something wrong with NSS',nss
 629       goto 13
 630    12 print *,'something wrong with header'
 631
 632    13 NCON=ICON-1
 633
 634 #endif
 635       call flush(iout)
 636       jj_old=1
 637       open (icbase,file=bprotfiles,status="unknown",&
 638          form="unformatted",access="direct",recl=lenrec)
 639 ! Read conformations from binary DA files (one per batch) and write them to
 640 ! a binary DA scratchfile.
 641       jj=0
 642       jjj=0
 643 #ifdef MPI
 644       write (liczba,'(bz,i3.3)') me
 645       IF (ME.EQ.MASTER) THEN
 646 ! Only the master reads the database; it'll send it to the other procs
 647 ! through a ring.
 648 #endif
 649         t_acq = tcpu()
 650         icount=0
 651
 652         if (from_bx) then
 653
 654           open (intin,file=intinname,status="old",form="unformatted",&
 655                   access="direct",recl=lenrec_in)
 656
 657         else if (from_cx) then
 658 #if (defined(AIX) && !defined(JUBL))
 659           call xdrfopen_(ixdrf,intinname, "r", iret)
 660 #else
 661           call xdrfopen(ixdrf,intinname, "r", iret)
 662 #endif
 663           prec=10000.0
 664           write (iout,*) "xdrfopen: iret",iret
 665           if (iret.eq.0) then
 666             write (iout,*) "Error: coordinate file ",&
 667              intinname(:ilen(intinname))," does not exist."
 668             call flush(iout)
 669 #ifdef MPI
 670             call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
 671 #endif
 672             stop
 673           endif
 674         else
 675           write (iout,*) "Error: coordinate format not specified"
 676           call flush(iout)
 677 #ifdef MPI
 678           call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
 679 #else
 680           stop
 681 #endif
 682         endif
 683
 684 !#define DEBUG
 685 #ifdef DEBUG
 686         write (iout,*) "Opening file ",intinname(:ilen(intinname))
 687         write (iout,*) "lenrec",lenrec_in
 688         call flush(iout)
 689 #endif
 690 !#undef DEBUG
 691 !        write (iout,*) "maxconf",maxconf
 692         i=0
 693         do while (.true.)
 694            i=i+1
 695 !el           if (i.gt.maxconf) then
 696 !el             write (iout,*) "Error: too many conformations ",&
 697 !el              "(",maxconf,") maximum."
 698 !#ifdef MPI
 699 !el             call MPI_Abort(MPI_COMM_WORLD,errcode,ierror)
 700 !#endif
 701 !el             stop
 702 !el           endif
 703 !          write (iout,*) "i",i
 704 !          call flush(iout)
 705           if (from_bx) then
 706             read(intin,err=101,end=101) &
 707              ((csingle(l,k),l=1,3),k=1,nres),&
 708              ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 709              nss,(ihpb(k),jhpb(k),k=1,nss),&
 710              energy(jj+1),&
 711              entfac(jj+1),rmstb(jj+1),iscor
 712              do j=1,2*nres
 713                do k=1,3
 714                  c(k,j)=csingle(k,j)
 715                enddo
 716              enddo
 717           else
 718 #if (defined(AIX) && !defined(JUBL))
 719             call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
 720             if (iret.eq.0) goto 101
 721             call xdrfint_(ixdrf, nss, iret)
 722             if (iret.eq.0) goto 101
 723             do j=1,nss
 724               call xdrfint_(ixdrf, ihpb(j), iret)
 725               if (iret.eq.0) goto 101
 726               call xdrfint_(ixdrf, jhpb(j), iret)
 727               if (iret.eq.0) goto 101
 728             enddo
 729             call xdrffloat_(ixdrf,reini,iret)
 730             if (iret.eq.0) goto 101
 731             call xdrffloat_(ixdrf,refree,iret)
 732             if (iret.eq.0) goto 101
 733             call xdrffloat_(ixdrf,rmsdev,iret)
 734             if (iret.eq.0) goto 101
 735             call xdrfint_(ixdrf,iscor,iret)
 736             if (iret.eq.0) goto 101
 737 #else
 738 !            write (iout,*) "calling xdrf3dfcoord"
 739             call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
 740 !            write (iout,*) "iret",iret
 741 !            call flush(iout)
 742             if (iret.eq.0) goto 101
 743             call xdrfint(ixdrf, nss, iret)
 744 !            write (iout,*) "iret",iret
 745 !            write (iout,*) "nss",nss
 746             call flush(iout)
 747             if (iret.eq.0) goto 101
 748             do k=1,nss
 749               call xdrfint(ixdrf, ihpb(k), iret)
 750               if (iret.eq.0) goto 101
 751               call xdrfint(ixdrf, jhpb(k), iret)
 752               if (iret.eq.0) goto 101
 753             enddo
 754             call xdrffloat(ixdrf,reini,iret)
 755             if (iret.eq.0) goto 101
 756             call xdrffloat(ixdrf,refree,iret)
 757             if (iret.eq.0) goto 101
 758             call xdrffloat(ixdrf,rmsdev,iret)
 759             if (iret.eq.0) goto 101
 760             call xdrfint(ixdrf,iscor,iret)
 761             if (iret.eq.0) goto 101
 762 #endif
 763             energy(jj+1)=reini
 764             entfac(jj+1)=refree
 765             rmstb(jj+1)=rmsdev
 766             do k=1,nres
 767               do l=1,3
 768                 c(l,k)=csingle(l,k)
 769               enddo
 770             enddo
 771             do k=nnt,nct
 772               do l=1,3
 773                 c(l,nres+k)=csingle(l,nres+k-nnt+1)
 774               enddo
 775             enddo
 776           endif
 777 #ifdef DEBUG
 778           write (iout,'(5hREAD ,i5,3f15.4,i10)') &
 779            jj+1,energy(jj+1),entfac(jj+1),&
 780            rmstb(jj+1),iscor
 781           write (iout,*) "Conformation",jjj+1,jj+1
 782           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 783           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 784           call flush(iout)
 785 #endif
 786           call add_new_cconf(jjj,jj,jj_old,icount,Next)
 787         enddo
 788   101   continue
 789         write (iout,*) i-1," conformations read from DA file ",&
 790           intinname(:ilen(intinname))
 791         write (iout,*) jj," conformations read so far"
 792         if (from_bx) then
 793           close(intin)
 794         else
 795 #if (defined(AIX) && !defined(JUBL))
 796           call xdrfclose_(ixdrf, iret)
 797 #else
 798           call xdrfclose(ixdrf, iret)
 799 #endif
 800         endif
 801 #ifdef MPI
 802 !#ifdef DEBUG
 803         write (iout,*) "jj_old",jj_old," jj",jj
 804 !#endif
 805         call write_and_send_cconf(icount,jj_old,jj,Next)
 806         call MPI_Send(0,1,MPI_INTEGER,Next,570,&
 807                    MPI_COMM_WORLD,IERROR)
 808         jj_old=jj+1
 809 #else
 810         call write_and_send_cconf(icount,jj_old,jj,Next)
 811 #endif
 812         t_acq = tcpu() - t_acq
 813 #ifdef MPI
 814         write (iout,*) "Processor",me,&
 815           " time for conformation read/send",t_acq
 816       ELSE
 817 ! A worker gets the confs from the master and sends them to its neighbor
 818         t_acq = tcpu()
 819         call receive_and_pass_cconf(icount,jj_old,jj,&
 820           Previous,Next)
 821         t_acq = tcpu() - t_acq
 822       ENDIF
 823 #endif
 824       ncon=jj
 825 !      close(icbase)
 826       close(intin)
 827
 828       write(iout,*)"A total of",ncon," conformations read."
 829
 830       allocate(enetb(1:max_ene,ncon)) !(max_ene,maxstr_proc)
 831 #ifdef MPI
 832 ! Check if everyone has the same number of conformations
 833       call MPI_Allgather(ncon,1,MPI_INTEGER,&
 834         ntot_all(0),1,MPI_INTEGER,MPI_Comm_World,IERROR)
 835       lerr=.false.
 836       do i=0,nprocs-1
 837         if (i.ne.me) then
 838             if (ncon.ne.ntot_all(i)) then
 839               write (iout,*) "Number of conformations at processor",i,&
 840                " differs from that at processor",me,&
 841                ncon,ntot_all(i)
 842               lerr = .true.
 843             endif
 844         endif
 845       enddo
 846       if (lerr) then
 847         write (iout,*)
 848         write (iout,*) "Number of conformations read by processors"
 849         write (iout,*)
 850         do i=0,nprocs-1
 851           write (iout,'(8i10)') i,ntot_all(i)
 852         enddo
 853         write (iout,*) "Calculation terminated."
 854         call flush(iout)
 855         return 1
 856       endif
 857       return
 858 #endif
 859  1111 write(iout,*) "Error opening coordinate file ",&
 860        intinname(:ilen(intinname))
 861       call flush(iout)
 862       return 1
 863       end subroutine read_coords
 864 !------------------------------------------------------------------------------
 865       subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
 866
 867       use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
 868       use energy_data, only: itel,itype,dsc,max_ene,molnum
 869       use control_data, only: symetr
 870       use geometry, only: int_from_cart1
 871 !      implicit none
 872 !      include "DIMENSIONS"
 873 !      include "sizesclu.dat"
 874 !      include "COMMON.CLUSTER"
 875 !      include "COMMON.CONTROL"
 876 !      include "COMMON.CHAIN"
 877 !      include "COMMON.INTERACT"
 878 !      include "COMMON.LOCAL"
 879 !      include "COMMON.IOUNITS"
 880 !      include "COMMON.NAMES"
 881 !      include "COMMON.VAR"
 882 !      include "COMMON.SBRIDGE"
 883 !      include "COMMON.GEO"
 884       integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
 885         nn,nn1,inan,Next,itj,chalen,mnum
 886       real(kind=8) :: etot,energia(0:max_ene)
 887       jjj=jjj+1
 888       chalen=int((nct-nnt+2)/symetr)
 889       call int_from_cart1(.false.)
 890       do j=nnt+1,nct
 891         mnum=molnum(j)
 892         if (mnum.eq.5) cycle
 893         if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
 894         if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle
 895
 896         if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
 897          if (j.gt.2) then
 898           if (itel(j).ne.0 .and. itel(j-1).ne.0) then
 899           write (iout,*) "Conformation",jjj,jj+1
 900           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
 901        chalen
 902           write (iout,*) "The Cartesian geometry is:"
 903           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 904           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 905           write (iout,*) "The internal geometry is:"
 906           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 907           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 908           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 909           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 910           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 911           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 912           write (iout,*) &
 913             "This conformation WILL NOT be added to the database."
 914           return
 915           endif
 916          endif
 917         endif
 918       enddo
 919       do j=nnt,nct
 920         mnum=molnum(j)
 921         itj=itype(j,mnum)
 922         if (mnum.eq.5) cycle
 923         if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) &
 924             .and.  (vbld(nres+j)-dsc(iabs(itj))) &
 925                                   .gt.2.0d0) then
 926           write (iout,*) "Conformation",jjj,jj+1
 927           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
 928           write (iout,*) "The Cartesian geometry is:"
 929           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 930           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 931           write (iout,*) "The internal geometry is:"
 932           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 933           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 934           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 935           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 936           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 937           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 938           write (iout,*) &
 939             "This conformation WILL NOT be added to the database."
 940           return
 941         endif
 942       enddo
 943       do j=3,nres
 944         if (theta(j).le.0.0d0) then
 945           write (iout,*) &
 946             "Zero theta angle(s) in conformation",jjj,jj+1
 947           write (iout,*) "The Cartesian geometry is:"
 948           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 949           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 950           write (iout,*) "The internal geometry is:"
 951           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 952           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 953           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 954           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 955           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 956           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 957           write (iout,*) &
 958             "This conformation WILL NOT be added to the database."
 959           return
 960         endif
 961         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 962       enddo
 963       jj=jj+1
 964 #ifdef DEBUG
 965       write (iout,*) "Conformation",jjj,jj
 966       write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 967       write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 968       write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 969       write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 970       write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 971       write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 972       write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
 973       write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 974       write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 975       write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 976       write (iout,'(e15.5,16i5)') entfac(icount+1)
 977 !     &        iscore(icount+1,0)
 978 #endif
 979       icount=icount+1
 980       call store_cconf_from_file(jj,icount)
 981       if (icount.eq.maxstr_proc) then
 982 #ifdef DEBUG
 983         write (iout,* ) "jj_old",jj_old," jj",jj
 984 #endif
 985         call write_and_send_cconf(icount,jj_old,jj,Next)
 986         jj_old=jj+1
 987         icount=0
 988       endif
 989       return
 990       end subroutine add_new_cconf
 991 !------------------------------------------------------------------------------
 992       subroutine store_cconf_from_file(jj,icount)
 993
 994       use energy_data, only: ihpb,jhpb
 995 !      implicit none
 996 !      include "DIMENSIONS"
 997 !      include "sizesclu.dat"
 998 !      include "COMMON.CLUSTER"
 999 !      include "COMMON.CHAIN"
1000 !      include "COMMON.SBRIDGE"
1001 !      include "COMMON.INTERACT"
1002 !      include "COMMON.IOUNITS"
1003 !      include "COMMON.VAR"
1004       integer :: i,j,jj,icount
1005 ! Store the conformation that has been read in
1006       do i=1,2*nres
1007         do j=1,3
1008           allcart(j,i,icount)=c(j,i)
1009         enddo
1010       enddo
1011       nss_all(icount)=nss
1012       do i=1,nss
1013         ihpb_all(i,icount)=ihpb(i)
1014         jhpb_all(i,icount)=jhpb(i)
1015       enddo
1016       return
1017       end subroutine store_cconf_from_file
1018 !------------------------------------------------------------------------------
1019       subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1020
1021 !      implicit none
1022 !      include "DIMENSIONS"
1023 !      include "sizesclu.dat"
1024 #ifdef MPI
1025       use MPI_data
1026       include "mpif.h"
1027       integer :: IERROR
1028 !      include "COMMON.MPI"
1029 #endif
1030 !      include "COMMON.CHAIN"
1031 !      include "COMMON.SBRIDGE"
1032 !      include "COMMON.INTERACT"
1033 !      include "COMMON.IOUNITS"
1034 !      include "COMMON.CLUSTER"
1035 !      include "COMMON.VAR"
1036       integer :: icount,jj_old,jj,Next
1037 ! Write the structures to a scratch file
1038 #ifdef MPI
1039 ! Master sends the portion of conformations that have been read in to the neighbor
1040 #ifdef DEBUG
1041       write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1042       call flush(iout)
1043 #endif
1044       call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1045       call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1046           Next,571,MPI_COMM_WORLD,IERROR)
1047       call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1048           Next,572,MPI_COMM_WORLD,IERROR)
1049       call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1050           Next,573,MPI_COMM_WORLD,IERROR)
1051       call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1052           Next,577,MPI_COMM_WORLD,IERROR)
1053       call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1054           Next,579,MPI_COMM_WORLD,IERROR)
1055       call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1056           MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1057 #endif
1058       call dawrite_ccoords(jj_old,jj,icbase)
1059       return
1060       end subroutine write_and_send_cconf
1061 !------------------------------------------------------------------------------
1062 #ifdef MPI
1063       subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1064
1065       use MPI_data
1066 !      implicit none
1067 !      include "DIMENSIONS"
1068 !      include "sizesclu.dat"
1069       include "mpif.h"
1070       integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1071 !      include "COMMON.MPI"
1072 !      include "COMMON.CHAIN"
1073 !      include "COMMON.SBRIDGE"
1074 !      include "COMMON.INTERACT"
1075 !      include "COMMON.IOUNITS"
1076 !      include "COMMON.VAR"
1077 !      include "COMMON.GEO"
1078 !      include "COMMON.CLUSTER"
1079       integer :: i,j,k,icount,jj_old,jj,Previous,Next
1080       icount=1
1081 #ifdef DEBUG
1082       write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1083       call flush(iout)
1084 #endif
1085       do while (icount.gt.0)
1086       call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1087            STATUS,IERROR)
1088       call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1089            IERROR)
1090 #ifdef DEBUG
1091       write (iout,*) "Processor",me," icount",icount
1092 #endif
1093       if (icount.eq.0) return
1094       call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1095           Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1096       call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1097         Next,571,MPI_COMM_WORLD,IERROR)
1098       call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1099           Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1100       call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1101         Next,572,MPI_COMM_WORLD,IERROR)
1102       call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1103           Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1104       call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1105         Next,573,MPI_COMM_WORLD,IERROR)
1106       call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1107         Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1108       call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1109         Next,577,MPI_COMM_WORLD,IERROR)
1110       call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1111         Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1112       call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1113         Next,579,MPI_COMM_WORLD,IERROR)
1114       call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1115         MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1116       call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1117         MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1118       jj=jj_old+icount-1
1119       call dawrite_ccoords(jj_old,jj,icbase)
1120       jj_old=jj+1
1121 #ifdef DEBUG
1122       write (iout,*) "Processor",me," received",icount," conformations"
1123       do i=1,icount
1124         write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1125         write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1126         write (iout,'(e15.5,16i5)') entfac(i)
1127       enddo
1128 #endif
1129       enddo
1130       return
1131       end subroutine receive_and_pass_cconf
1132 #endif
1133 !------------------------------------------------------------------------------
1134       subroutine daread_ccoords(istart_conf,iend_conf)
1135
1136 !      implicit none
1137 !      include "DIMENSIONS"
1138 !      include "sizesclu.dat"
1139 #ifdef MPI
1140       use MPI_data
1141       include "mpif.h"
1142 !      include "COMMON.MPI"
1143 #endif
1144 !      include "COMMON.CHAIN"
1145 !      include "COMMON.CLUSTER"
1146 !      include "COMMON.IOUNITS"
1147 !      include "COMMON.INTERACT"
1148 !      include "COMMON.VAR"
1149 !      include "COMMON.SBRIDGE"
1150 !      include "COMMON.GEO"
1151       integer :: istart_conf,iend_conf
1152       integer :: i,j,ij,ii,iii
1153       integer :: len
1154       character(len=16) :: form,acc
1155       character(len=32) :: nam
1156 !
1157 ! Read conformations off a DA scratchfile.
1158 !
1159 #ifdef DEBUG
1160       write (iout,*) "DAREAD_COORDS"
1161       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1162       inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1163       write (iout,*) "len=",len," form=",form," acc=",acc
1164       write (iout,*) "nam=",nam
1165       call flush(iout)
1166 #endif
1167       do ii=istart_conf,iend_conf
1168         ij = ii - istart_conf + 1
1169         iii=list_conf(ii)
1170 #ifdef DEBUG
1171         write (iout,*) "Reading binary file, record",iii," ii",ii
1172         call flush(iout)
1173 #endif
1174         read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1175           ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1176           nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1177           entfac(ii),rmstb(ii)
1178 #ifdef DEBUG
1179         write (iout,*) ii,iii,ij,entfac(ii)
1180         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1181         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1182           i=nnt+nres,nct+nres)
1183         write (iout,'(2e15.5)') entfac(ij)
1184         write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1185           jhpb_all(i,ij),i=1,nss)
1186         call flush(iout)
1187 #endif
1188       enddo
1189       return
1190       end subroutine daread_ccoords
1191 !------------------------------------------------------------------------------
1192       subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1193
1194 !      implicit none
1195 !      include "DIMENSIONS"
1196 !      include "sizesclu.dat"
1197 #ifdef MPI
1198       use MPI_data
1199       include "mpif.h"
1200 !      include "COMMON.MPI"
1201 #endif
1202 !      include "COMMON.CHAIN"
1203 !      include "COMMON.INTERACT"
1204 !      include "COMMON.IOUNITS"
1205 !      include "COMMON.VAR"
1206 !      include "COMMON.SBRIDGE"
1207 !      include "COMMON.GEO"
1208 !      include "COMMON.CLUSTER"
1209       integer :: istart_conf,iend_conf
1210       integer :: i,j,ii,ij,iii,unit_out
1211       integer :: len
1212       character(len=16) :: form,acc
1213       character(len=32) :: nam
1214 !
1215 ! Write conformations to a DA scratchfile.
1216 !
1217 #ifdef DEBUG
1218       write (iout,*) "DAWRITE_COORDS"
1219       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1220       write (iout,*) "lenrec",lenrec
1221       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1222       write (iout,*) "len=",len," form=",form," acc=",acc
1223       write (iout,*) "nam=",nam
1224       call flush(iout)
1225 #endif
1226       do ii=istart_conf,iend_conf
1227         iii=list_conf(ii)
1228         ij = ii - istart_conf + 1
1229 #ifdef DEBUG
1230         write (iout,*) "Writing binary file, record",iii," ii",ii
1231         call flush(iout)
1232 #endif
1233         write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1234           ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1235           nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1236           entfac(ii),rmstb(ii)
1237 #ifdef DEBUG
1238         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1239         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1240          nct+nres)
1241         write (iout,'(2e15.5)') entfac(ij)
1242         write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1243          nss_all(ij))
1244         call flush(iout)
1245 #endif
1246       enddo
1247       return
1248       end subroutine dawrite_ccoords
1249 !-----------------------------------------------------------------------------
1250 ! readrtns.F
1251 !-----------------------------------------------------------------------------
1252       subroutine read_control
1253 !
1254 ! Read molecular data
1255 !
1256       use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
1257       use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1258                  iscode,symetr,punch_dist,print_dist,nstart,nend,&
1259                  caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
1260                  r_cut_ele
1261 !      implicit none
1262 !      include 'DIMENSIONS'
1263 !      include 'sizesclu.dat'
1264 !      include 'COMMON.IOUNITS'
1265 !      include 'COMMON.TIME1'
1266 !      include 'COMMON.SBRIDGE'
1267 !      include 'COMMON.CONTROL'
1268 !      include 'COMMON.CLUSTER'
1269 !      include 'COMMON.CHAIN'
1270 !      include 'COMMON.HEADER'
1271 !      include 'COMMON.FFIELD'
1272 !      include 'COMMON.FREE'
1273 !      include 'COMMON.INTERACT'
1274       character(len=320) :: controlcard !,ucase
1275 !#ifdef MPL
1276 !      include 'COMMON.INFO'
1277 !#endif
1278       integer :: i
1279
1280       read (INP,'(a80)') titel
1281       call card_concat(controlcard,.true.)
1282
1283       call readi(controlcard,'NRES',nres_molec(1),0)
1284       call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
1285       call readi(controlcard,"NRES_CAT",nres_molec(5),0)
1286       nres=0
1287       do i=1,5
1288        nres=nres_molec(i)+nres
1289       enddo
1290
1291 !      call alloc_clust_arrays
1292       allocate(rcutoff(max_cut+1)) !(max_cut+1)
1293       allocate(beta_h(maxT)) !(maxT)
1294       allocate(mult(nres)) !(maxres)
1295
1296
1297       call readi(controlcard,'RESCALE',rescale_mode,2)
1298       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1299       write (iout,*) "DISTCHAINMAX",distchainmax
1300       call readi(controlcard,'PDBOUT',outpdb,0)
1301       call readi(controlcard,'MOL2OUT',outmol2,0)
1302       refstr=(index(controlcard,'REFSTR').gt.0)
1303       write (iout,*) "REFSTR",refstr
1304       pdbref=(index(controlcard,'PDBREF').gt.0)
1305       iscode=index(controlcard,'ONE_LETTER')
1306       tree=(index(controlcard,'MAKE_TREE').gt.0)
1307       min_var=(index(controlcard,'MINVAR').gt.0)
1308       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1309       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1310       call readi(controlcard,'NCUT',ncut,1)
1311       call readi(controlcard,'SYM',symetr,1)
1312       write (iout,*) 'sym', symetr
1313       call readi(controlcard,'NSTART',nstart,0)
1314       call reada(controlcard,'BOXX',boxxsize,100.0d0)
1315       call reada(controlcard,'BOXY',boxysize,100.0d0)
1316       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1317 !      ions=index(controlcard,"IONS").gt.0
1318       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
1319       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1320       write(iout,*) "R_CUT_ELE=",r_cut_ele
1321       call readi(controlcard,'NEND',nend,0)
1322       call reada(controlcard,'ECUT',ecut,10.0d0)
1323       call reada(controlcard,'PROB',prob_limit,0.99d0)
1324       write (iout,*) "Probability limit",prob_limit
1325       lgrp=(index(controlcard,'LGRP').gt.0)
1326       caonly=(index(controlcard,'CA_ONLY').gt.0)
1327       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1328       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1329       call readi(controlcard,'IOPT',iopt,2)
1330       lside = index(controlcard,"SIDE").gt.0
1331       efree = index(controlcard,"EFREE").gt.0
1332       call readi(controlcard,'NTEMP',nT,1)
1333       write (iout,*) "nT",nT
1334 !el      call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1335       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1336       write (iout,*) "nT",nT
1337       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1338       do i=1,nT
1339         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1340       enddo
1341       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1342       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1343       lprint_int=index(controlcard,"PRINT_INT") .gt.0
1344       if (min_var) iopt=1
1345       return
1346       end subroutine read_control
1347 !-----------------------------------------------------------------------------
1348       subroutine molread
1349 !
1350 ! Read molecular data.
1351 !
1352       use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
1353       use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1354 !                 wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1355 !                 wturn3,wturn4,wturn6,wvdwpp,weights
1356       use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1357                  indpdb
1358       use geometry, only: chainbuild,alloc_geo_arrays
1359       use energy, only: alloc_ener_arrays
1360       use io_base, only: rescode
1361       use control, only: setup_var,init_int_table
1362       use conform_compar, only: contact
1363 !      implicit none
1364 !      include 'DIMENSIONS'
1365 !      include 'COMMON.IOUNITS'
1366 !      include 'COMMON.GEO'
1367 !      include 'COMMON.VAR'
1368 !      include 'COMMON.INTERACT'
1369 !      include 'COMMON.LOCAL'
1370 !      include 'COMMON.NAMES'
1371 !      include 'COMMON.CHAIN'
1372 !      include 'COMMON.FFIELD'
1373 !      include 'COMMON.SBRIDGE'
1374 !      include 'COMMON.HEADER'
1375 !      include 'COMMON.CONTROL'
1376 !      include 'COMMON.CONTACTS'
1377 !      include 'COMMON.TIME1'
1378 !#ifdef MPL
1379 !      include 'COMMON.INFO'
1380 !#endif
1381       character(len=4) ::  sequence(nres) !(maxres)
1382       character(len=800) :: weightcard
1383 !      integer rescode
1384       real(kind=8) :: x(6*nres) !(maxvar)
1385       integer  :: itype_pdb(nres) !(maxres)
1386 !      logical seq_comp
1387       integer :: i,j,kkk,mnum
1388 !
1389 ! Body
1390 !
1391 !el      allocate(weights(n_ene))
1392       allocate(weights(max_ene))
1393       call alloc_geo_arrays
1394       call alloc_ener_arrays
1395 !-----------------------------
1396       allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1397       allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1398       allocate(itype(nres+2,5)) !(maxres)
1399       allocate(molnum(nres+1))
1400       allocate(itel(nres+2))
1401
1402       do i=1,2*nres+2
1403         do j=1,3
1404           c(j,i)=0
1405           dc(j,i)=0
1406         enddo
1407       enddo
1408       do j=1,5
1409       do i=1,nres+2
1410         itype(i,j)=0
1411         itel(i)=0
1412       enddo
1413       enddo
1414 !--------------------------
1415 ! Read weights of the subsequent energy terms.
1416       call card_concat(weightcard,.true.)
1417       call reada(weightcard,'WSC',wsc,1.0d0)
1418       call reada(weightcard,'WLONG',wsc,wsc)
1419       call reada(weightcard,'WSCP',wscp,1.0d0)
1420       call reada(weightcard,'WELEC',welec,1.0D0)
1421       call reada(weightcard,'WVDWPP',wvdwpp,welec)
1422       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1423       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1424       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1425       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1426       call reada(weightcard,'WTURN3',wturn3,1.0D0)
1427       call reada(weightcard,'WTURN4',wturn4,1.0D0)
1428       call reada(weightcard,'WTURN6',wturn6,1.0D0)
1429       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1430       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1431       call reada(weightcard,'WBOND',wbond,1.0D0)
1432       call reada(weightcard,'WTOR',wtor,1.0D0)
1433       call reada(weightcard,'WTORD',wtor_d,1.0D0)
1434       call reada(weightcard,'WANG',wang,1.0D0)
1435       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1436       call reada(weightcard,'SCAL14',scal14,0.4D0)
1437       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1438       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1439       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1440       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
1441       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
1442       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
1443       if (index(weightcard,'SOFT').gt.0) ipot=6
1444 ! 12/1/95 Added weight for the multi-body term WCORR
1445       call reada(weightcard,'WCORRH',wcorr,1.0D0)
1446       if (wcorr4.gt.0.0d0) wcorr=wcorr4
1447       weights(1)=wsc
1448       weights(2)=wscp
1449       weights(3)=welec
1450       weights(4)=wcorr
1451       weights(5)=wcorr5
1452       weights(6)=wcorr6
1453       weights(7)=wel_loc
1454       weights(8)=wturn3
1455       weights(9)=wturn4
1456       weights(10)=wturn6
1457       weights(11)=wang
1458       weights(12)=wscloc
1459       weights(13)=wtor
1460       weights(14)=wtor_d
1461       weights(15)=wstrain
1462       weights(16)=wvdwpp
1463       weights(17)=wbond
1464       weights(18)=scal14
1465 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
1466       weights(21)=wsccor
1467       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1468         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1469         wturn4,wturn6,wsccor
1470    10 format (/'Energy-term weights (unscaled):'// &
1471        'WSCC=   ',f10.6,' (SC-SC)'/ &
1472        'WSCP=   ',f10.6,' (SC-p)'/ &
1473        'WELEC=  ',f10.6,' (p-p electr)'/ &
1474        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1475        'WBOND=  ',f10.6,' (stretching)'/ &
1476        'WANG=   ',f10.6,' (bending)'/ &
1477        'WSCLOC= ',f10.6,' (SC local)'/ &
1478        'WTOR=   ',f10.6,' (torsional)'/ &
1479        'WTORD=  ',f10.6,' (double torsional)'/ &
1480        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1481        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1482        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1483        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1484        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1485        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1486        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1487        'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1488        'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1489
1490       if (wcorr4.gt.0.0d0) then
1491         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1492          'between contact pairs of peptide groups'
1493         write (iout,'(2(a,f5.3/))') &
1494         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1495         'Range of quenching the correlation terms:',2*delt_corr
1496       else if (wcorr.gt.0.0d0) then
1497         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1498          'between contact pairs of peptide groups'
1499       endif
1500       write (iout,'(a,f8.3)') &
1501         'Scaling factor of 1,4 SC-p interactions:',scal14
1502       write (iout,'(a,f8.3)') &
1503         'General scaling factor of SC-p interactions:',scalscp
1504       r0_corr=cutoff_corr-delt_corr
1505       do i=1,20
1506         aad(i,1)=scalscp*aad(i,1)
1507         aad(i,2)=scalscp*aad(i,2)
1508         bad(i,1)=scalscp*bad(i,1)
1509         bad(i,2)=scalscp*bad(i,2)
1510       enddo
1511
1512       call flush(iout)
1513       print *,'indpdb=',indpdb,' pdbref=',pdbref
1514
1515 ! Read sequence if not taken from the pdb file.
1516       if (iscode.gt.0) then
1517         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1518       else
1519         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1520       endif
1521 ! Convert sequence to numeric code
1522       do i=1,nres_molec(1)
1523         mnum=1
1524         molnum(i)=1
1525         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1526       enddo
1527       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
1528         mnum=2
1529         molnum(i)=2
1530         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1531       enddo
1532       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
1533         mnum=5
1534         molnum(i)=5
1535         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1536       enddo
1537       print *,nres
1538       print '(20i4)',(itype(i,molnum(i)),i=1,nres)
1539
1540       do i=1,nres-1
1541 #ifdef PROCOR
1542         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
1543             .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
1544 #else
1545         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
1546 #endif
1547           itel(i)=0
1548 #ifdef PROCOR
1549         else if (iabs(itype(i+1,1)).ne.20) then
1550 #else
1551         else if (iabs(itype(i,1)).ne.20) then
1552 #endif
1553           itel(i)=1
1554         else
1555           itel(i)=2
1556         endif
1557       enddo
1558       write (iout,*) "ITEL"
1559       do i=1,nres-1
1560         write (iout,*) i,itype(i,molnum(i)),itel(i)
1561       enddo
1562
1563       print *,'Call Read_Bridge.'
1564       call read_bridge
1565       nnt=1
1566       nct=nres
1567       print *,'NNT=',NNT,' NCT=',NCT
1568       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1569       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1570       if (nstart.lt.nnt) nstart=nnt
1571       if (nend.gt.nct .or. nend.eq.0) nend=nct
1572       write (iout,*) "nstart",nstart," nend",nend
1573       nres0=nres
1574 !      if (pdbref) then
1575 !        read(inp,'(a)') pdbfile
1576 !        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1577 !        open(ipdbin,file=pdbfile,status='old',err=33)
1578 !        goto 34
1579 !  33    write (iout,'(a)') 'Error opening PDB file.'
1580 !        stop
1581 !  34    continue
1582 !        print *,'Begin reading pdb data'
1583 !        call readpdb
1584 !        print *,'Finished reading pdb data'
1585 !        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1586 !        do i=1,nres
1587 !          itype_pdb(i)=itype(i)
1588 !        enddo
1589 !        close (ipdbin)
1590 !        write (iout,'(a,i3)') 'nsup=',nsup
1591 !        nstart_seq=nnt
1592 !        if (nsup.le.(nct-nnt+1)) then
1593 !          do i=0,nct-nnt+1-nsup
1594 !            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1595 !              nstart_seq=nnt+i
1596 !              goto 111
1597 !            endif
1598 !          enddo
1599 !          write (iout,'(a)')
1600 !     &            'Error - sequences to be superposed do not match.'
1601 !          stop
1602 !        else
1603 !          do i=0,nsup-(nct-nnt+1)
1604 !            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1605 !     &      then
1606 !              nstart_sup=nstart_sup+i
1607 !              nsup=nct-nnt+1
1608 !              goto 111
1609 !            endif
1610 !          enddo
1611 !          write (iout,'(a)')
1612 !     &            'Error - sequences to be superposed do not match.'
1613 !        endif
1614 !  111   continue
1615 !        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1616 !     &                 ' nstart_seq=',nstart_seq
1617 !      endif
1618 write(iout,*)"przed ini_int_tab"
1619       call init_int_table
1620 write(iout,*)"po ini_int_tab"
1621 write(iout,*)"przed setup var"
1622       call setup_var
1623 write(iout,*)"po setup var"
1624       write (iout,*) "molread: REFSTR",refstr
1625       if (refstr) then
1626         if (.not.pdbref) then
1627           call read_angles(inp,*38)
1628           goto 39
1629    38     write (iout,'(a)') 'Error reading reference structure.'
1630 #ifdef MPL
1631           call mp_stopall(Error_Msg)
1632 #else
1633           stop 'Error reading reference structure'
1634 #endif
1635    39     call chainbuild
1636           nstart_sup=nnt
1637           nstart_seq=nnt
1638           nsup=nct-nnt+1
1639           kkk=1
1640           do i=1,2*nres
1641             do j=1,3
1642               cref(j,i,kkk)=c(j,i)
1643             enddo
1644           enddo
1645         endif
1646         call contact(.true.,ncont_ref,icont_ref)
1647       endif
1648       return
1649       end subroutine molread
1650 !-----------------------------------------------------------------------------
1651       subroutine openunits
1652 !      implicit none
1653 !      include 'DIMENSIONS'
1654       use control_data, only: from_cx,from_bx,from_cart
1655 #ifdef MPI
1656       use MPI_data
1657       include "mpif.h"
1658       character(len=3) :: liczba
1659 !      include "COMMON.MPI"
1660 #endif
1661 !      include 'COMMON.IOUNITS'
1662 !      include 'COMMON.CONTROL'
1663       integer :: lenpre,lenpot !,ilen
1664 !      external ilen
1665       character(len=16) :: cformat,cprint
1666 !      character(len=16) ucase
1667       integer :: lenint,lenout
1668       call getenv('INPUT',prefix)
1669       call getenv('OUTPUT',prefout)
1670       call getenv('INTIN',prefintin)
1671       call getenv('COORD',cformat)
1672       call getenv('PRINTCOOR',cprint)
1673       call getenv('SCRATCHDIR',scratchdir)
1674       from_bx=.true.
1675       from_cx=.false.
1676       if (index(ucase(cformat),'CX').gt.0) then
1677         from_cx=.true.
1678         from_bx=.false.
1679       endif
1680       from_cart=.true.
1681       lenpre=ilen(prefix)
1682       lenout=ilen(prefout)
1683       lenint=ilen(prefintin)
1684 ! Get the names and open the input files
1685       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1686 #ifdef MPI
1687       write (liczba,'(bz,i3.3)') me
1688       outname=prefout(:lenout)//'_clust.out_'//liczba
1689 #else
1690       outname=prefout(:lenout)//'_clust.out'
1691 #endif
1692       if (from_bx) then
1693         intinname=prefintin(:lenint)//'.bx'
1694       else if (from_cx) then
1695         intinname=prefintin(:lenint)//'.cx'
1696       else
1697         intinname=prefintin(:lenint)//'.int'
1698       endif
1699       rmsname=prefintin(:lenint)//'.rms'
1700       open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1701              status='unknown')
1702       open (jrms,file=rmsname,status='unknown')
1703       open(iout,file=outname,status='unknown')
1704 ! Get parameter filenames and open the parameter files.
1705       call getenv('BONDPAR',bondname)
1706       open (ibond,file=bondname,status='old')
1707       call getenv('THETPAR',thetname)
1708       open (ithep,file=thetname,status='old')
1709       call getenv('ROTPAR',rotname)
1710       open (irotam,file=rotname,status='old')
1711       call getenv('TORPAR',torname)
1712       open (itorp,file=torname,status='old')
1713       call getenv('TORDPAR',tordname)
1714       open (itordp,file=tordname,status='old')
1715       call getenv('FOURIER',fouriername)
1716       open (ifourier,file=fouriername,status='old')
1717       call getenv('ELEPAR',elename)
1718       open (ielep,file=elename,status='old')
1719       call getenv('SIDEPAR',sidename)
1720       open (isidep,file=sidename,status='old')
1721       call getenv('SIDEP',sidepname)
1722       open (isidep1,file=sidepname,status="old")
1723       call getenv('SCCORPAR',sccorname)
1724       open (isccor,file=sccorname,status="old")
1725       call getenv('THETPAR_NUCL',thetname_nucl)
1726       open (ithep_nucl,file=thetname_nucl,status='old')
1727       call getenv('ROTPAR_NUCL',rotname_nucl)
1728       open (irotam_nucl,file=rotname_nucl,status='old')
1729       call getenv('TORPAR_NUCL',torname_nucl)
1730       open (itorp_nucl,file=torname_nucl,status='old')
1731       call getenv('TORDPAR_NUCL',tordname_nucl)
1732       open (itordp_nucl,file=tordname_nucl,status='old')
1733       call getenv('SIDEPAR_NUCL',sidename_nucl)
1734       open (isidep_nucl,file=sidename_nucl,status='old')
1735       call getenv('SIDEPAR_SCBASE',sidename_scbase)
1736       open (isidep_scbase,file=sidename_scbase,status='old')
1737       call getenv('PEPPAR_PEPBASE',pepname_pepbase)
1738       open (isidep_pepbase,file=pepname_pepbase,status='old')
1739       call getenv('SCPAR_PHOSPH',pepname_scpho)
1740       open (isidep_scpho,file=pepname_scpho,status='old')
1741       call getenv('PEPPAR_PHOSPH',pepname_peppho)
1742       open (isidep_peppho,file=pepname_peppho,status='old')
1743
1744
1745       call getenv('LIPTRANPAR',liptranname)
1746       open (iliptranpar,file=liptranname,status='old')
1747       call getenv('TUBEPAR',tubename)
1748       open (itube,file=tubename,status='old')
1749       call getenv('IONPAR',ionname)
1750       open (iion,file=ionname,status='old')
1751
1752 #ifndef OLDSCP
1753 !
1754 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1755 ! Use -DOLDSCP to use hard-coded constants instead.
1756 !
1757       call getenv('SCPPAR',scpname)
1758       open (iscpp,file=scpname,status='old')
1759 #endif
1760       return
1761       end subroutine openunits
1762 !-----------------------------------------------------------------------------
1763 ! geomout.F
1764 !-----------------------------------------------------------------------------
1765       subroutine pdboutC(etot,rmsd,tytul)
1766
1767       use energy_data, only: ihpb,jhpb,itype,molnum
1768 !      implicit real*8 (a-h,o-z)
1769 !      include 'DIMENSIONS'
1770 !      include 'COMMON.CONTROL'
1771 !      include 'COMMON.CHAIN'
1772 !      include 'COMMON.INTERACT'
1773 !      include 'COMMON.NAMES'
1774 !      include 'COMMON.IOUNITS'
1775 !      include 'COMMON.HEADER'
1776 !      include 'COMMON.SBRIDGE'
1777 !      include 'COMMON.TEMPFAC'
1778       character(len=50) :: tytul
1779       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1780                                                'G','H','I','J'/)
1781       integer :: ica(nres)
1782       real(kind=8) :: etot,rmsd
1783       integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
1784       real(kind=8) :: boxxxx(3)
1785       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1786         ' ENERGY ',etot,' RMS ',rmsd
1787       iatom=0
1788       ichain=1
1789       ires=0
1790       boxxxshift(1)=int(c(1,nnt)/boxxsize)
1791 !      if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
1792       boxxxshift(2)=int(c(2,nnt)/boxzsize)
1793 !      if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
1794       boxxxshift(3)=int(c(3,nnt)/boxzsize)
1795 !      if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
1796
1797       boxxxx(1)=boxxsize
1798       boxxxx(2)=boxysize
1799       boxxxx(3)=boxzsize
1800
1801       do i=nnt,nct
1802         mnum=molnum(i)
1803         iti=itype(i,mnum)
1804         if (iti.eq.ntyp1_molec(mnum)) then
1805           ichain=ichain+1
1806           ires=0
1807           write (ipdb,'(a)') 'TER'
1808         else
1809         ires=ires+1
1810         iatom=iatom+1
1811         ica(i)=iatom
1812         if (mnum.eq.5) then
1813            do j=1,3
1814             if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then
1815            c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j)
1816             else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then
1817            c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j)
1818              else
1819            c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j)
1820             endif
1821            enddo
1822         endif
1823         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
1824            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
1825         if ((iti.ne.10).and.(mnum.ne.5)) then
1826           iatom=iatom+1
1827           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
1828             ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
1829         endif
1830         endif
1831       enddo
1832       write (ipdb,'(a)') 'TER'
1833       do i=nnt,nct-1
1834         mnum=molnum(i)
1835         mnum1=molnum(i+1)
1836         if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
1837         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1838           write (ipdb,30) ica(i),ica(i+1)
1839         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1840           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
1841         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
1842           write (ipdb,30) ica(i),ica(i)+1
1843         endif
1844       enddo
1845       if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
1846         write (ipdb,30) ica(nct),ica(nct)+1
1847       endif
1848       do i=1,nss
1849         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1850       enddo
1851       write (ipdb,'(a6)') 'ENDMDL'
1852   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1853   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1854   30  FORMAT ('CONECT',8I5)
1855       return
1856       end subroutine pdboutC
1857 !-----------------------------------------------------------------------------
1858       subroutine cartout(igr,i,etot,free,rmsd,plik)
1859 !     implicit real*8 (a-h,o-z)
1860 !     include 'DIMENSIONS'
1861 !     include 'sizesclu.dat'
1862 !     include 'COMMON.IOUNITS'
1863 !     include 'COMMON.CHAIN'
1864 !     include 'COMMON.VAR'
1865 !     include 'COMMON.LOCAL'
1866 !     include 'COMMON.INTERACT'
1867 !     include 'COMMON.NAMES'
1868 !     include 'COMMON.GEO'
1869 !     include 'COMMON.CLUSTER'
1870       integer :: igr,i,j,k
1871       real(kind=8) :: etot,free,rmsd
1872       character(len=80) :: plik
1873       open (igeom,file=plik,position='append')
1874       write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
1875       write (igeom,'(i4,$)') &
1876         nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
1877       write (igeom,'(i10)') iscore(i)
1878       write (igeom,'(8f10.5)') &
1879         ((allcart(k,j,i),k=1,3),j=1,nres),&
1880         ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
1881       return
1882       end subroutine cartout
1883 !------------------------------------------------------------------------------
1884 !      subroutine alloc_clust_arrays(n_conf)
1885
1886 !      integer :: n_conf
1887 !COMMON.CLUSTER
1888 !      common /clu/
1889 !      allocate(diss(maxdist)) !(maxdist)
1890 !el      allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
1891 !      allocatable :: enetb !(max_ene,maxstr_proc)
1892 !el      allocate(entfac(maxconf)) !(maxconf)
1893 !      allocatable :: totfree_gr !(maxgr)
1894 !el      allocate(rcutoff(max_cut+1)) !(max_cut+1)
1895 !      common /clu1/
1896 !      allocatable :: licz,iass !(maxgr)
1897 !      allocatable :: nconf !(maxgr,maxingr)
1898 !      allocatable :: iass_tot !(maxgr,max_cut)
1899 !      allocatable :: list_conf !(maxconf)
1900 !      common /alles/
1901 !el      allocatable :: allcart !(3,maxres2,maxstr_proc)
1902 !el      allocate(rmstb(maxconf)) !(maxconf)
1903 !el      allocate(mult(nres)) !(maxres)
1904 !el      allocatable :: nss_all !(maxstr_proc)
1905 !el      allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
1906 !      allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
1907 !COMMON.TEMPFAC
1908 !      common /factemp/
1909 !      allocatable :: tempfac !(2,maxres)
1910 !COMMON.FREE
1911 !      common /free/
1912 !el      allocate(beta_h(maxT)) !(maxT)
1913
1914 !      end subroutine alloc_clust_arrays
1915 !-----------------------------------------------------------------------------
1916 !-----------------------------------------------------------------------------
1917       end module io_clust