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
   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
 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
 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         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
 892          if (j.gt.2) then
 893           if (itel(j).ne.0 .and. itel(j-1).ne.0) then
 894           write (iout,*) "Conformation",jjj,jj+1
 895           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
 896        chalen
 897           write (iout,*) "The Cartesian geometry is:"
 898           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 899           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 900           write (iout,*) "The internal geometry is:"
 901           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 902           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 903           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 904           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 905           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 906           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 907           write (iout,*) &
 908             "This conformation WILL NOT be added to the database."
 909           return
 910           endif
 911          endif
 912         endif
 913       enddo
 914       do j=nnt,nct
 915         itj=itype(j)
 916         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))) &
 917                                   .gt.2.0d0) then
 918           write (iout,*) "Conformation",jjj,jj+1
 919           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
 920           write (iout,*) "The Cartesian geometry is:"
 921           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 922           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 923           write (iout,*) "The internal geometry is:"
 924           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 925           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 926           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 927           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 928           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 929           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 930           write (iout,*) &
 931             "This conformation WILL NOT be added to the database."
 932           return
 933         endif
 934       enddo
 935       do j=3,nres
 936         if (theta(j).le.0.0d0) then
 937           write (iout,*) &
 938             "Zero theta angle(s) in conformation",jjj,jj+1
 939           write (iout,*) "The Cartesian geometry is:"
 940           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 941           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 942           write (iout,*) "The internal geometry is:"
 943           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 944           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 945           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 946           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 947           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 948           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 949           write (iout,*) &
 950             "This conformation WILL NOT be added to the database."
 951           return
 952         endif
 953         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 954       enddo
 955       jj=jj+1
 956 #ifdef DEBUG
 957       write (iout,*) "Conformation",jjj,jj
 958       write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 959       write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 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)') (vbld(k+nres),k=nnt,nct)
 965       write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 966       write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 967       write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 968       write (iout,'(e15.5,16i5)') entfac(icount+1)
 969 !     &        iscore(icount+1,0)
 970 #endif
 971       icount=icount+1
 972       call store_cconf_from_file(jj,icount)
 973       if (icount.eq.maxstr_proc) then
 974 #ifdef DEBUG
 975         write (iout,* ) "jj_old",jj_old," jj",jj
 976 #endif
 977         call write_and_send_cconf(icount,jj_old,jj,Next)
 978         jj_old=jj+1
 979         icount=0
 980       endif
 981       return
 982       end subroutine add_new_cconf
 983 !------------------------------------------------------------------------------
 984       subroutine store_cconf_from_file(jj,icount)
 985
 986       use energy_data, only: ihpb,jhpb
 987 !      implicit none
 988 !      include "DIMENSIONS"
 989 !      include "sizesclu.dat"
 990 !      include "COMMON.CLUSTER"
 991 !      include "COMMON.CHAIN"
 992 !      include "COMMON.SBRIDGE"
 993 !      include "COMMON.INTERACT"
 994 !      include "COMMON.IOUNITS"
 995 !      include "COMMON.VAR"
 996       integer :: i,j,jj,icount
 997 ! Store the conformation that has been read in
 998       do i=1,2*nres
 999         do j=1,3
1000           allcart(j,i,icount)=c(j,i)
1001         enddo
1002       enddo
1003       nss_all(icount)=nss
1004       do i=1,nss
1005         ihpb_all(i,icount)=ihpb(i)
1006         jhpb_all(i,icount)=jhpb(i)
1007       enddo
1008       return
1009       end subroutine store_cconf_from_file
1010 !------------------------------------------------------------------------------
1011       subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1012
1013 !      implicit none
1014 !      include "DIMENSIONS"
1015 !      include "sizesclu.dat"
1016 #ifdef MPI
1017       use MPI_data
1018       include "mpif.h"
1019       integer :: IERROR
1020 !      include "COMMON.MPI"
1021 #endif
1022 !      include "COMMON.CHAIN"
1023 !      include "COMMON.SBRIDGE"
1024 !      include "COMMON.INTERACT"
1025 !      include "COMMON.IOUNITS"
1026 !      include "COMMON.CLUSTER"
1027 !      include "COMMON.VAR"
1028       integer :: icount,jj_old,jj,Next
1029 ! Write the structures to a scratch file
1030 #ifdef MPI
1031 ! Master sends the portion of conformations that have been read in to the neighbor
1032 #ifdef DEBUG
1033       write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1034       call flush(iout)
1035 #endif
1036       call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1037       call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1038           Next,571,MPI_COMM_WORLD,IERROR)
1039       call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1040           Next,572,MPI_COMM_WORLD,IERROR)
1041       call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1042           Next,573,MPI_COMM_WORLD,IERROR)
1043       call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1044           Next,577,MPI_COMM_WORLD,IERROR)
1045       call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1046           Next,579,MPI_COMM_WORLD,IERROR)
1047       call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1048           MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1049 #endif
1050       call dawrite_ccoords(jj_old,jj,icbase)
1051       return
1052       end subroutine write_and_send_cconf
1053 !------------------------------------------------------------------------------
1054 #ifdef MPI
1055       subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1056
1057       use MPI_data
1058 !      implicit none
1059 !      include "DIMENSIONS"
1060 !      include "sizesclu.dat"
1061       include "mpif.h"
1062       integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1063 !      include "COMMON.MPI"
1064 !      include "COMMON.CHAIN"
1065 !      include "COMMON.SBRIDGE"
1066 !      include "COMMON.INTERACT"
1067 !      include "COMMON.IOUNITS"
1068 !      include "COMMON.VAR"
1069 !      include "COMMON.GEO"
1070 !      include "COMMON.CLUSTER"
1071       integer :: i,j,k,icount,jj_old,jj,Previous,Next
1072       icount=1
1073 #ifdef DEBUG
1074       write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1075       call flush(iout)
1076 #endif
1077       do while (icount.gt.0)
1078       call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1079            STATUS,IERROR)
1080       call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1081            IERROR)
1082 #ifdef DEBUG
1083       write (iout,*) "Processor",me," icount",icount
1084 #endif
1085       if (icount.eq.0) return
1086       call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1087           Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1088       call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1089         Next,571,MPI_COMM_WORLD,IERROR)
1090       call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1091           Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1092       call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1093         Next,572,MPI_COMM_WORLD,IERROR)
1094       call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1095           Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1096       call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1097         Next,573,MPI_COMM_WORLD,IERROR)
1098       call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1099         Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1100       call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1101         Next,577,MPI_COMM_WORLD,IERROR)
1102       call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1103         Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1104       call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1105         Next,579,MPI_COMM_WORLD,IERROR)
1106       call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1107         MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1108       call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1109         MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1110       jj=jj_old+icount-1
1111       call dawrite_ccoords(jj_old,jj,icbase)
1112       jj_old=jj+1
1113 #ifdef DEBUG
1114       write (iout,*) "Processor",me," received",icount," conformations"
1115       do i=1,icount
1116         write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1117         write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1118         write (iout,'(e15.5,16i5)') entfac(i)
1119       enddo
1120 #endif
1121       enddo
1122       return
1123       end subroutine receive_and_pass_cconf
1124 #endif
1125 !------------------------------------------------------------------------------
1126       subroutine daread_ccoords(istart_conf,iend_conf)
1127
1128 !      implicit none
1129 !      include "DIMENSIONS"
1130 !      include "sizesclu.dat"
1131 #ifdef MPI
1132       use MPI_data
1133       include "mpif.h"
1134 !      include "COMMON.MPI"
1135 #endif
1136 !      include "COMMON.CHAIN"
1137 !      include "COMMON.CLUSTER"
1138 !      include "COMMON.IOUNITS"
1139 !      include "COMMON.INTERACT"
1140 !      include "COMMON.VAR"
1141 !      include "COMMON.SBRIDGE"
1142 !      include "COMMON.GEO"
1143       integer :: istart_conf,iend_conf
1144       integer :: i,j,ij,ii,iii
1145       integer :: len
1146       character(len=16) :: form,acc
1147       character(len=32) :: nam
1148 !
1149 ! Read conformations off a DA scratchfile.
1150 !
1151 #ifdef DEBUG
1152       write (iout,*) "DAREAD_COORDS"
1153       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1154       inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1155       write (iout,*) "len=",len," form=",form," acc=",acc
1156       write (iout,*) "nam=",nam
1157       call flush(iout)
1158 #endif
1159       do ii=istart_conf,iend_conf
1160         ij = ii - istart_conf + 1
1161         iii=list_conf(ii)
1162 #ifdef DEBUG
1163         write (iout,*) "Reading binary file, record",iii," ii",ii
1164         call flush(iout)
1165 #endif
1166         read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1167           ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1168           nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1169           entfac(ii),rmstb(ii)
1170 #ifdef DEBUG
1171         write (iout,*) ii,iii,ij,entfac(ii)
1172         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1173         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1174           i=nnt+nres,nct+nres)
1175         write (iout,'(2e15.5)') entfac(ij)
1176         write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1177           jhpb_all(i,ij),i=1,nss)
1178         call flush(iout)
1179 #endif
1180       enddo
1181       return
1182       end subroutine daread_ccoords
1183 !------------------------------------------------------------------------------
1184       subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1185
1186 !      implicit none
1187 !      include "DIMENSIONS"
1188 !      include "sizesclu.dat"
1189 #ifdef MPI
1190       use MPI_data
1191       include "mpif.h"
1192 !      include "COMMON.MPI"
1193 #endif
1194 !      include "COMMON.CHAIN"
1195 !      include "COMMON.INTERACT"
1196 !      include "COMMON.IOUNITS"
1197 !      include "COMMON.VAR"
1198 !      include "COMMON.SBRIDGE"
1199 !      include "COMMON.GEO"
1200 !      include "COMMON.CLUSTER"
1201       integer :: istart_conf,iend_conf
1202       integer :: i,j,ii,ij,iii,unit_out
1203       integer :: len
1204       character(len=16) :: form,acc
1205       character(len=32) :: nam
1206 !
1207 ! Write conformations to a DA scratchfile.
1208 !
1209 #ifdef DEBUG
1210       write (iout,*) "DAWRITE_COORDS"
1211       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1212       write (iout,*) "lenrec",lenrec
1213       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1214       write (iout,*) "len=",len," form=",form," acc=",acc
1215       write (iout,*) "nam=",nam
1216       call flush(iout)
1217 #endif
1218       do ii=istart_conf,iend_conf
1219         iii=list_conf(ii)
1220         ij = ii - istart_conf + 1
1221 #ifdef DEBUG
1222         write (iout,*) "Writing binary file, record",iii," ii",ii
1223         call flush(iout)
1224 #endif
1225         write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1226           ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1227           nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1228           entfac(ii),rmstb(ii)
1229 #ifdef DEBUG
1230         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1231         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1232          nct+nres)
1233         write (iout,'(2e15.5)') entfac(ij)
1234         write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1235          nss_all(ij))
1236         call flush(iout)
1237 #endif
1238       enddo
1239       return
1240       end subroutine dawrite_ccoords
1241 !-----------------------------------------------------------------------------
1242 ! readrtns.F
1243 !-----------------------------------------------------------------------------
1244       subroutine read_control
1245 !
1246 ! Read molecular data
1247 !
1248       use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
1249       use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1250                  iscode,symetr,punch_dist,print_dist,nstart,nend,&
1251                  caonly,iopt,efree,lprint_cart,lprint_int
1252 !      implicit none
1253 !      include 'DIMENSIONS'
1254 !      include 'sizesclu.dat'
1255 !      include 'COMMON.IOUNITS'
1256 !      include 'COMMON.TIME1'
1257 !      include 'COMMON.SBRIDGE'
1258 !      include 'COMMON.CONTROL'
1259 !      include 'COMMON.CLUSTER'
1260 !      include 'COMMON.CHAIN'
1261 !      include 'COMMON.HEADER'
1262 !      include 'COMMON.FFIELD'
1263 !      include 'COMMON.FREE'
1264 !      include 'COMMON.INTERACT'
1265       character(len=320) :: controlcard !,ucase
1266 !#ifdef MPL
1267 !      include 'COMMON.INFO'
1268 !#endif
1269       integer :: i
1270
1271       read (INP,'(a80)') titel
1272       call card_concat(controlcard,.true.)
1273
1274       call readi(controlcard,'NRES',nres,0)
1275
1276 !      call alloc_clust_arrays
1277       allocate(rcutoff(max_cut+1)) !(max_cut+1)
1278       allocate(beta_h(maxT)) !(maxT)
1279       allocate(mult(nres)) !(maxres)
1280
1281
1282       call readi(controlcard,'RESCALE',rescale_mode,2)
1283       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1284       write (iout,*) "DISTCHAINMAX",distchainmax
1285       call readi(controlcard,'PDBOUT',outpdb,0)
1286       call readi(controlcard,'MOL2OUT',outmol2,0)
1287       refstr=(index(controlcard,'REFSTR').gt.0)
1288       write (iout,*) "REFSTR",refstr
1289       pdbref=(index(controlcard,'PDBREF').gt.0)
1290       iscode=index(controlcard,'ONE_LETTER')
1291       tree=(index(controlcard,'MAKE_TREE').gt.0)
1292       min_var=(index(controlcard,'MINVAR').gt.0)
1293       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1294       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1295       call readi(controlcard,'NCUT',ncut,1)
1296       call readi(controlcard,'SYM',symetr,1)
1297       write (iout,*) 'sym', symetr
1298       call readi(controlcard,'NSTART',nstart,0)
1299       call readi(controlcard,'NEND',nend,0)
1300       call reada(controlcard,'ECUT',ecut,10.0d0)
1301       call reada(controlcard,'PROB',prob_limit,0.99d0)
1302       write (iout,*) "Probability limit",prob_limit
1303       lgrp=(index(controlcard,'LGRP').gt.0)
1304       caonly=(index(controlcard,'CA_ONLY').gt.0)
1305       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1306       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1307       call readi(controlcard,'IOPT',iopt,2)
1308       lside = index(controlcard,"SIDE").gt.0
1309       efree = index(controlcard,"EFREE").gt.0
1310       call readi(controlcard,'NTEMP',nT,1)
1311       write (iout,*) "nT",nT
1312 !el      call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1313       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1314       write (iout,*) "nT",nT
1315       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1316       do i=1,nT
1317         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1318       enddo
1319       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1320       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1321       lprint_int=index(controlcard,"PRINT_INT") .gt.0
1322       if (min_var) iopt=1
1323       return
1324       end subroutine read_control
1325 !-----------------------------------------------------------------------------
1326       subroutine molread
1327 !
1328 ! Read molecular data.
1329 !
1330       use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
1331       use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1332 !                 wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1333 !                 wturn3,wturn4,wturn6,wvdwpp,weights
1334       use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1335                  indpdb
1336       use geometry, only: chainbuild,alloc_geo_arrays
1337       use energy, only: alloc_ener_arrays
1338       use control, only: rescode,setup_var,init_int_table
1339       use conform_compar, only: contact
1340 !      implicit none
1341 !      include 'DIMENSIONS'
1342 !      include 'COMMON.IOUNITS'
1343 !      include 'COMMON.GEO'
1344 !      include 'COMMON.VAR'
1345 !      include 'COMMON.INTERACT'
1346 !      include 'COMMON.LOCAL'
1347 !      include 'COMMON.NAMES'
1348 !      include 'COMMON.CHAIN'
1349 !      include 'COMMON.FFIELD'
1350 !      include 'COMMON.SBRIDGE'
1351 !      include 'COMMON.HEADER'
1352 !      include 'COMMON.CONTROL'
1353 !      include 'COMMON.CONTACTS'
1354 !      include 'COMMON.TIME1'
1355 !#ifdef MPL
1356 !      include 'COMMON.INFO'
1357 !#endif
1358       character(len=4) ::  sequence(nres) !(maxres)
1359       character(len=800) :: weightcard
1360 !      integer rescode
1361       real(kind=8) :: x(6*nres) !(maxvar)
1362       integer  :: itype_pdb(nres) !(maxres)
1363 !      logical seq_comp
1364       integer :: i,j,kkk
1365 !
1366 ! Body
1367 !
1368 !el      allocate(weights(n_ene))
1369       allocate(weights(max_ene))
1370       call alloc_geo_arrays
1371       call alloc_ener_arrays
1372 !-----------------------------
1373       allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1374       allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1375       allocate(itype(nres+2)) !(maxres)
1376       allocate(itel(nres+2))
1377
1378       do i=1,2*nres+2
1379         do j=1,3
1380           c(j,i)=0
1381           dc(j,i)=0
1382         enddo
1383       enddo
1384       do i=1,nres+2
1385         itype(i)=0
1386         itel(i)=0
1387       enddo
1388 !--------------------------
1389 ! Read weights of the subsequent energy terms.
1390       call card_concat(weightcard,.true.)
1391       call reada(weightcard,'WSC',wsc,1.0d0)
1392       call reada(weightcard,'WLONG',wsc,wsc)
1393       call reada(weightcard,'WSCP',wscp,1.0d0)
1394       call reada(weightcard,'WELEC',welec,1.0D0)
1395       call reada(weightcard,'WVDWPP',wvdwpp,welec)
1396       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1397       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1398       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1399       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1400       call reada(weightcard,'WTURN3',wturn3,1.0D0)
1401       call reada(weightcard,'WTURN4',wturn4,1.0D0)
1402       call reada(weightcard,'WTURN6',wturn6,1.0D0)
1403       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1404       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1405       call reada(weightcard,'WBOND',wbond,1.0D0)
1406       call reada(weightcard,'WTOR',wtor,1.0D0)
1407       call reada(weightcard,'WTORD',wtor_d,1.0D0)
1408       call reada(weightcard,'WANG',wang,1.0D0)
1409       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1410       call reada(weightcard,'SCAL14',scal14,0.4D0)
1411       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1412       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1413       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1414       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
1415       if (index(weightcard,'SOFT').gt.0) ipot=6
1416 ! 12/1/95 Added weight for the multi-body term WCORR
1417       call reada(weightcard,'WCORRH',wcorr,1.0D0)
1418       if (wcorr4.gt.0.0d0) wcorr=wcorr4
1419       weights(1)=wsc
1420       weights(2)=wscp
1421       weights(3)=welec
1422       weights(4)=wcorr
1423       weights(5)=wcorr5
1424       weights(6)=wcorr6
1425       weights(7)=wel_loc
1426       weights(8)=wturn3
1427       weights(9)=wturn4
1428       weights(10)=wturn6
1429       weights(11)=wang
1430       weights(12)=wscloc
1431       weights(13)=wtor
1432       weights(14)=wtor_d
1433       weights(15)=wstrain
1434       weights(16)=wvdwpp
1435       weights(17)=wbond
1436       weights(18)=scal14
1437 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
1438       weights(21)=wsccor
1439       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1440         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1441         wturn4,wturn6,wsccor
1442    10 format (/'Energy-term weights (unscaled):'// &
1443        'WSCC=   ',f10.6,' (SC-SC)'/ &
1444        'WSCP=   ',f10.6,' (SC-p)'/ &
1445        'WELEC=  ',f10.6,' (p-p electr)'/ &
1446        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1447        'WBOND=  ',f10.6,' (stretching)'/ &
1448        'WANG=   ',f10.6,' (bending)'/ &
1449        'WSCLOC= ',f10.6,' (SC local)'/ &
1450        'WTOR=   ',f10.6,' (torsional)'/ &
1451        'WTORD=  ',f10.6,' (double torsional)'/ &
1452        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1453        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1454        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1455        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1456        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1457        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1458        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1459        'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1460        'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1461
1462       if (wcorr4.gt.0.0d0) then
1463         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1464          'between contact pairs of peptide groups'
1465         write (iout,'(2(a,f5.3/))') &
1466         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1467         'Range of quenching the correlation terms:',2*delt_corr
1468       else if (wcorr.gt.0.0d0) then
1469         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1470          'between contact pairs of peptide groups'
1471       endif
1472       write (iout,'(a,f8.3)') &
1473         'Scaling factor of 1,4 SC-p interactions:',scal14
1474       write (iout,'(a,f8.3)') &
1475         'General scaling factor of SC-p interactions:',scalscp
1476       r0_corr=cutoff_corr-delt_corr
1477       do i=1,20
1478         aad(i,1)=scalscp*aad(i,1)
1479         aad(i,2)=scalscp*aad(i,2)
1480         bad(i,1)=scalscp*bad(i,1)
1481         bad(i,2)=scalscp*bad(i,2)
1482       enddo
1483
1484       call flush(iout)
1485       print *,'indpdb=',indpdb,' pdbref=',pdbref
1486
1487 ! Read sequence if not taken from the pdb file.
1488       if (iscode.gt.0) then
1489         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1490       else
1491         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1492       endif
1493 ! Convert sequence to numeric code
1494       do i=1,nres
1495         itype(i)=rescode(i,sequence(i),iscode)
1496       enddo
1497       print *,nres
1498       print '(20i4)',(itype(i),i=1,nres)
1499
1500       do i=1,nres
1501 #ifdef PROCOR
1502         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
1503 #else
1504         if (itype(i).eq.ntyp1) then
1505 #endif
1506           itel(i)=0
1507 #ifdef PROCOR
1508         else if (iabs(itype(i+1)).ne.20) then
1509 #else
1510         else if (iabs(itype(i)).ne.20) then
1511 #endif
1512           itel(i)=1
1513         else
1514           itel(i)=2
1515         endif
1516       enddo
1517       write (iout,*) "ITEL"
1518       do i=1,nres-1
1519         write (iout,*) i,itype(i),itel(i)
1520       enddo
1521
1522       print *,'Call Read_Bridge.'
1523       call read_bridge
1524       nnt=1
1525       nct=nres
1526       print *,'NNT=',NNT,' NCT=',NCT
1527       if (itype(1).eq.ntyp1) nnt=2
1528       if (itype(nres).eq.ntyp1) nct=nct-1
1529       if (nstart.lt.nnt) nstart=nnt
1530       if (nend.gt.nct .or. nend.eq.0) nend=nct
1531       write (iout,*) "nstart",nstart," nend",nend
1532       nres0=nres
1533 !      if (pdbref) then
1534 !        read(inp,'(a)') pdbfile
1535 !        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1536 !        open(ipdbin,file=pdbfile,status='old',err=33)
1537 !        goto 34
1538 !  33    write (iout,'(a)') 'Error opening PDB file.'
1539 !        stop
1540 !  34    continue
1541 !        print *,'Begin reading pdb data'
1542 !        call readpdb
1543 !        print *,'Finished reading pdb data'
1544 !        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1545 !        do i=1,nres
1546 !          itype_pdb(i)=itype(i)
1547 !        enddo
1548 !        close (ipdbin)
1549 !        write (iout,'(a,i3)') 'nsup=',nsup
1550 !        nstart_seq=nnt
1551 !        if (nsup.le.(nct-nnt+1)) then
1552 !          do i=0,nct-nnt+1-nsup
1553 !            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1554 !              nstart_seq=nnt+i
1555 !              goto 111
1556 !            endif
1557 !          enddo
1558 !          write (iout,'(a)')
1559 !     &            'Error - sequences to be superposed do not match.'
1560 !          stop
1561 !        else
1562 !          do i=0,nsup-(nct-nnt+1)
1563 !            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1564 !     &      then
1565 !              nstart_sup=nstart_sup+i
1566 !              nsup=nct-nnt+1
1567 !              goto 111
1568 !            endif
1569 !          enddo
1570 !          write (iout,'(a)')
1571 !     &            'Error - sequences to be superposed do not match.'
1572 !        endif
1573 !  111   continue
1574 !        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1575 !     &                 ' nstart_seq=',nstart_seq
1576 !      endif
1577 write(iout,*)"przed ini_int_tab"
1578       call init_int_table
1579 write(iout,*)"po ini_int_tab"
1580 write(iout,*)"przed setup var"
1581       call setup_var
1582 write(iout,*)"po setup var"
1583       write (iout,*) "molread: REFSTR",refstr
1584       if (refstr) then
1585         if (.not.pdbref) then
1586           call read_angles(inp,*38)
1587           goto 39
1588    38     write (iout,'(a)') 'Error reading reference structure.'
1589 #ifdef MPL
1590           call mp_stopall(Error_Msg)
1591 #else
1592           stop 'Error reading reference structure'
1593 #endif
1594    39     call chainbuild
1595           nstart_sup=nnt
1596           nstart_seq=nnt
1597           nsup=nct-nnt+1
1598           kkk=1
1599           do i=1,2*nres
1600             do j=1,3
1601               cref(j,i,kkk)=c(j,i)
1602             enddo
1603           enddo
1604         endif
1605         call contact(.true.,ncont_ref,icont_ref)
1606       endif
1607       return
1608       end subroutine molread
1609 !-----------------------------------------------------------------------------
1610       subroutine openunits
1611 !      implicit none
1612 !      include 'DIMENSIONS'
1613       use control_data, only: from_cx,from_bx,from_cart
1614 #ifdef MPI
1615       use MPI_data
1616       include "mpif.h"
1617       character(len=3) :: liczba
1618 !      include "COMMON.MPI"
1619 #endif
1620 !      include 'COMMON.IOUNITS'
1621 !      include 'COMMON.CONTROL'
1622       integer :: lenpre,lenpot !,ilen
1623 !      external ilen
1624       character(len=16) :: cformat,cprint
1625 !      character(len=16) ucase
1626       integer :: lenint,lenout
1627       call getenv('INPUT',prefix)
1628       call getenv('OUTPUT',prefout)
1629       call getenv('INTIN',prefintin)
1630       call getenv('COORD',cformat)
1631       call getenv('PRINTCOOR',cprint)
1632       call getenv('SCRATCHDIR',scratchdir)
1633       from_bx=.true.
1634       from_cx=.false.
1635       if (index(ucase(cformat),'CX').gt.0) then
1636         from_cx=.true.
1637         from_bx=.false.
1638       endif
1639       from_cart=.true.
1640       lenpre=ilen(prefix)
1641       lenout=ilen(prefout)
1642       lenint=ilen(prefintin)
1643 ! Get the names and open the input files
1644       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1645 #ifdef MPI
1646       write (liczba,'(bz,i3.3)') me
1647       outname=prefout(:lenout)//'_clust.out_'//liczba
1648 #else
1649       outname=prefout(:lenout)//'_clust.out'
1650 #endif
1651       if (from_bx) then
1652         intinname=prefintin(:lenint)//'.bx'
1653       else if (from_cx) then
1654         intinname=prefintin(:lenint)//'.cx'
1655       else
1656         intinname=prefintin(:lenint)//'.int'
1657       endif
1658       rmsname=prefintin(:lenint)//'.rms'
1659       open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1660              status='unknown')
1661       open (jrms,file=rmsname,status='unknown')
1662       open(iout,file=outname,status='unknown')
1663 ! Get parameter filenames and open the parameter files.
1664       call getenv('BONDPAR',bondname)
1665       open (ibond,file=bondname,status='old')
1666       call getenv('THETPAR',thetname)
1667       open (ithep,file=thetname,status='old')
1668       call getenv('ROTPAR',rotname)
1669       open (irotam,file=rotname,status='old')
1670       call getenv('TORPAR',torname)
1671       open (itorp,file=torname,status='old')
1672       call getenv('TORDPAR',tordname)
1673       open (itordp,file=tordname,status='old')
1674       call getenv('FOURIER',fouriername)
1675       open (ifourier,file=fouriername,status='old')
1676       call getenv('ELEPAR',elename)
1677       open (ielep,file=elename,status='old')
1678       call getenv('SIDEPAR',sidename)
1679       open (isidep,file=sidename,status='old')
1680       call getenv('SIDEP',sidepname)
1681       open (isidep1,file=sidepname,status="old")
1682       call getenv('SCCORPAR',sccorname)
1683       open (isccor,file=sccorname,status="old")
1684 #ifndef OLDSCP
1685 !
1686 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1687 ! Use -DOLDSCP to use hard-coded constants instead.
1688 !
1689       call getenv('SCPPAR',scpname)
1690       open (iscpp,file=scpname,status='old')
1691 #endif
1692       return
1693       end subroutine openunits
1694 !-----------------------------------------------------------------------------
1695 ! geomout.F
1696 !-----------------------------------------------------------------------------
1697       subroutine pdboutC(etot,rmsd,tytul)
1698
1699       use energy_data, only: ihpb,jhpb,itype
1700 !      implicit real*8 (a-h,o-z)
1701 !      include 'DIMENSIONS'
1702 !      include 'COMMON.CONTROL'
1703 !      include 'COMMON.CHAIN'
1704 !      include 'COMMON.INTERACT'
1705 !      include 'COMMON.NAMES'
1706 !      include 'COMMON.IOUNITS'
1707 !      include 'COMMON.HEADER'
1708 !      include 'COMMON.SBRIDGE'
1709 !      include 'COMMON.TEMPFAC'
1710       character(len=50) :: tytul
1711       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1712                                                'G','H','I','J'/)
1713       integer :: ica(nres)
1714       real(kind=8) :: etot,rmsd
1715       integer :: iatom,ichain,ires,i,j,iti
1716
1717       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1718         ' ENERGY ',etot,' RMS ',rmsd
1719       iatom=0
1720       ichain=1
1721       ires=0
1722       do i=nnt,nct
1723         iti=itype(i)
1724         if (iti.eq.ntyp1) then
1725           ichain=ichain+1
1726           ires=0
1727           write (ipdb,'(a)') 'TER'
1728         else
1729         ires=ires+1
1730         iatom=iatom+1
1731         ica(i)=iatom
1732         write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
1733            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
1734         if (iti.ne.10) then
1735           iatom=iatom+1
1736           write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
1737             ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
1738         endif
1739         endif
1740       enddo
1741       write (ipdb,'(a)') 'TER'
1742       do i=nnt,nct-1
1743         if (itype(i).eq.ntyp1) cycle
1744         if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
1745           write (ipdb,30) ica(i),ica(i+1)
1746         else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
1747           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
1748         else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
1749           write (ipdb,30) ica(i),ica(i)+1
1750         endif
1751       enddo
1752       if (itype(nct).ne.10) then
1753         write (ipdb,30) ica(nct),ica(nct)+1
1754       endif
1755       do i=1,nss
1756         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1757       enddo
1758       write (ipdb,'(a6)') 'ENDMDL'
1759   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1760   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1761   30  FORMAT ('CONECT',8I5)
1762       return
1763       end subroutine pdboutC
1764 !-----------------------------------------------------------------------------
1765       subroutine cartout(igr,i,etot,free,rmsd,plik)
1766 !     implicit real*8 (a-h,o-z)
1767 !     include 'DIMENSIONS'
1768 !     include 'sizesclu.dat'
1769 !     include 'COMMON.IOUNITS'
1770 !     include 'COMMON.CHAIN'
1771 !     include 'COMMON.VAR'
1772 !     include 'COMMON.LOCAL'
1773 !     include 'COMMON.INTERACT'
1774 !     include 'COMMON.NAMES'
1775 !     include 'COMMON.GEO'
1776 !     include 'COMMON.CLUSTER'
1777       integer :: igr,i,j,k
1778       real(kind=8) :: etot,free,rmsd
1779       character(len=80) :: plik
1780       open (igeom,file=plik,position='append')
1781       write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
1782       write (igeom,'(i4,$)') &
1783         nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
1784       write (igeom,'(i10)') iscore(i)
1785       write (igeom,'(8f10.5)') &
1786         ((allcart(k,j,i),k=1,3),j=1,nres),&
1787         ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
1788       return
1789       end subroutine cartout
1790 !------------------------------------------------------------------------------
1791 !      subroutine alloc_clust_arrays(n_conf)
1792
1793 !      integer :: n_conf
1794 !COMMON.CLUSTER
1795 !      common /clu/
1796 !      allocate(diss(maxdist)) !(maxdist)
1797 !el      allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
1798 !      allocatable :: enetb !(max_ene,maxstr_proc)
1799 !el      allocate(entfac(maxconf)) !(maxconf)
1800 !      allocatable :: totfree_gr !(maxgr)
1801 !el      allocate(rcutoff(max_cut+1)) !(max_cut+1)
1802 !      common /clu1/
1803 !      allocatable :: licz,iass !(maxgr)
1804 !      allocatable :: nconf !(maxgr,maxingr)
1805 !      allocatable :: iass_tot !(maxgr,max_cut)
1806 !      allocatable :: list_conf !(maxconf)
1807 !      common /alles/
1808 !el      allocatable :: allcart !(3,maxres2,maxstr_proc)
1809 !el      allocate(rmstb(maxconf)) !(maxconf)
1810 !el      allocate(mult(nres)) !(maxres)
1811 !el      allocatable :: nss_all !(maxstr_proc)
1812 !el      allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
1813 !      allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
1814 !COMMON.TEMPFAC
1815 !      common /factemp/
1816 !      allocatable :: tempfac !(2,maxres)
1817 !COMMON.FREE
1818 !      common /free/
1819 !el      allocate(beta_h(maxT)) !(maxT)
1820
1821 !      end subroutine alloc_clust_arrays
1822 !-----------------------------------------------------------------------------
1823 !-----------------------------------------------------------------------------
1824       end module io_clust