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