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.eq.5) 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.eq.5) 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 ! 12/1/95 Added weight for the multi-body term WCORR
1452       call reada(weightcard,'WCORRH',wcorr,1.0D0)
1453       if (wcorr4.gt.0.0d0) wcorr=wcorr4
1454       weights(1)=wsc
1455       weights(2)=wscp
1456       weights(3)=welec
1457       weights(4)=wcorr
1458       weights(5)=wcorr5
1459       weights(6)=wcorr6
1460       weights(7)=wel_loc
1461       weights(8)=wturn3
1462       weights(9)=wturn4
1463       weights(10)=wturn6
1464       weights(11)=wang
1465       weights(12)=wscloc
1466       weights(13)=wtor
1467       weights(14)=wtor_d
1468       weights(15)=wstrain
1469       weights(16)=wvdwpp
1470       weights(17)=wbond
1471       weights(18)=scal14
1472 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
1473       weights(21)=wsccor
1474       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1475         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1476         wturn4,wturn6,wsccor
1477    10 format (/'Energy-term weights (unscaled):'// &
1478        'WSCC=   ',f10.6,' (SC-SC)'/ &
1479        'WSCP=   ',f10.6,' (SC-p)'/ &
1480        'WELEC=  ',f10.6,' (p-p electr)'/ &
1481        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1482        'WBOND=  ',f10.6,' (stretching)'/ &
1483        'WANG=   ',f10.6,' (bending)'/ &
1484        'WSCLOC= ',f10.6,' (SC local)'/ &
1485        'WTOR=   ',f10.6,' (torsional)'/ &
1486        'WTORD=  ',f10.6,' (double torsional)'/ &
1487        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1488        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1489        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1490        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1491        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1492        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1493        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1494        'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1495        'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1496
1497       if (wcorr4.gt.0.0d0) then
1498         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1499          'between contact pairs of peptide groups'
1500         write (iout,'(2(a,f5.3/))') &
1501         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1502         'Range of quenching the correlation terms:',2*delt_corr
1503       else if (wcorr.gt.0.0d0) then
1504         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1505          'between contact pairs of peptide groups'
1506       endif
1507       write (iout,'(a,f8.3)') &
1508         'Scaling factor of 1,4 SC-p interactions:',scal14
1509       write (iout,'(a,f8.3)') &
1510         'General scaling factor of SC-p interactions:',scalscp
1511       r0_corr=cutoff_corr-delt_corr
1512       do i=1,20
1513         aad(i,1)=scalscp*aad(i,1)
1514         aad(i,2)=scalscp*aad(i,2)
1515         bad(i,1)=scalscp*bad(i,1)
1516         bad(i,2)=scalscp*bad(i,2)
1517       enddo
1518
1519       call flush(iout)
1520       print *,'indpdb=',indpdb,' pdbref=',pdbref
1521
1522 ! Read sequence if not taken from the pdb file.
1523       if (iscode.gt.0) then
1524         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1525       else
1526         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1527       endif
1528 ! Convert sequence to numeric code
1529       do i=1,nres_molec(1)
1530         mnum=1
1531         molnum(i)=1
1532         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1533       enddo
1534       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
1535         mnum=2
1536         molnum(i)=2
1537         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1538       enddo
1539       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
1540         mnum=5
1541         molnum(i)=5
1542         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1543       enddo
1544       print *,nres
1545       print '(20i4)',(itype(i,molnum(i)),i=1,nres)
1546
1547       do i=1,nres-1
1548 #ifdef PROCOR
1549         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
1550             .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
1551 #else
1552         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
1553 #endif
1554           itel(i)=0
1555 #ifdef PROCOR
1556         else if (iabs(itype(i+1,1)).ne.20) then
1557 #else
1558         else if (iabs(itype(i,1)).ne.20) then
1559 #endif
1560           itel(i)=1
1561         else
1562           itel(i)=2
1563         endif
1564       enddo
1565       write (iout,*) "ITEL"
1566       do i=1,nres-1
1567         write (iout,*) i,itype(i,molnum(i)),itel(i)
1568       enddo
1569
1570       print *,'Call Read_Bridge.'
1571       call read_bridge
1572       nnt=1
1573       nct=nres
1574       print *,'NNT=',NNT,' NCT=',NCT
1575       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1576       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1577       if (nstart.lt.nnt) nstart=nnt
1578       if (nend.gt.nct .or. nend.eq.0) nend=nct
1579       write (iout,*) "nstart",nstart," nend",nend
1580       nres0=nres
1581 !      if (pdbref) then
1582 !        read(inp,'(a)') pdbfile
1583 !        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1584 !        open(ipdbin,file=pdbfile,status='old',err=33)
1585 !        goto 34
1586 !  33    write (iout,'(a)') 'Error opening PDB file.'
1587 !        stop
1588 !  34    continue
1589 !        print *,'Begin reading pdb data'
1590 !        call readpdb
1591 !        print *,'Finished reading pdb data'
1592 !        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1593 !        do i=1,nres
1594 !          itype_pdb(i)=itype(i)
1595 !        enddo
1596 !        close (ipdbin)
1597 !        write (iout,'(a,i3)') 'nsup=',nsup
1598 !        nstart_seq=nnt
1599 !        if (nsup.le.(nct-nnt+1)) then
1600 !          do i=0,nct-nnt+1-nsup
1601 !            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1602 !              nstart_seq=nnt+i
1603 !              goto 111
1604 !            endif
1605 !          enddo
1606 !          write (iout,'(a)')
1607 !     &            'Error - sequences to be superposed do not match.'
1608 !          stop
1609 !        else
1610 !          do i=0,nsup-(nct-nnt+1)
1611 !            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1612 !     &      then
1613 !              nstart_sup=nstart_sup+i
1614 !              nsup=nct-nnt+1
1615 !              goto 111
1616 !            endif
1617 !          enddo
1618 !          write (iout,'(a)')
1619 !     &            'Error - sequences to be superposed do not match.'
1620 !        endif
1621 !  111   continue
1622 !        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1623 !     &                 ' nstart_seq=',nstart_seq
1624 !      endif
1625 write(iout,*)"przed ini_int_tab"
1626       call init_int_table
1627 write(iout,*)"po ini_int_tab"
1628 write(iout,*)"przed setup var"
1629       call setup_var
1630 write(iout,*)"po setup var"
1631       write (iout,*) "molread: REFSTR",refstr
1632       if (refstr) then
1633         if (.not.pdbref) then
1634           call read_angles(inp,*38)
1635           goto 39
1636    38     write (iout,'(a)') 'Error reading reference structure.'
1637 #ifdef MPL
1638           call mp_stopall(Error_Msg)
1639 #else
1640           stop 'Error reading reference structure'
1641 #endif
1642    39     call chainbuild
1643           nstart_sup=nnt
1644           nstart_seq=nnt
1645           nsup=nct-nnt+1
1646           kkk=1
1647           do i=1,2*nres
1648             do j=1,3
1649               cref(j,i,kkk)=c(j,i)
1650             enddo
1651           enddo
1652         endif
1653         call contact(.true.,ncont_ref,icont_ref)
1654       endif
1655       return
1656       end subroutine molread
1657 !-----------------------------------------------------------------------------
1658       subroutine openunits
1659 !      implicit none
1660 !      include 'DIMENSIONS'
1661       use control_data, only: from_cx,from_bx,from_cart
1662 #ifdef MPI
1663       use MPI_data
1664       include "mpif.h"
1665       character(len=3) :: liczba
1666 !      include "COMMON.MPI"
1667 #endif
1668 !      include 'COMMON.IOUNITS'
1669 !      include 'COMMON.CONTROL'
1670       integer :: lenpre,lenpot !,ilen
1671 !      external ilen
1672       character(len=16) :: cformat,cprint
1673 !      character(len=16) ucase
1674       integer :: lenint,lenout
1675       call getenv('INPUT',prefix)
1676       call getenv('OUTPUT',prefout)
1677       call getenv('INTIN',prefintin)
1678       call getenv('COORD',cformat)
1679       call getenv('PRINTCOOR',cprint)
1680       call getenv('SCRATCHDIR',scratchdir)
1681       from_bx=.true.
1682       from_cx=.false.
1683       if (index(ucase(cformat),'CX').gt.0) then
1684         from_cx=.true.
1685         from_bx=.false.
1686       endif
1687       from_cart=.true.
1688       lenpre=ilen(prefix)
1689       lenout=ilen(prefout)
1690       lenint=ilen(prefintin)
1691 ! Get the names and open the input files
1692       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1693 #ifdef MPI
1694       write (liczba,'(bz,i3.3)') me
1695       outname=prefout(:lenout)//'_clust.out_'//liczba
1696 #else
1697       outname=prefout(:lenout)//'_clust.out'
1698 #endif
1699       if (from_bx) then
1700         intinname=prefintin(:lenint)//'.bx'
1701       else if (from_cx) then
1702         intinname=prefintin(:lenint)//'.cx'
1703       else
1704         intinname=prefintin(:lenint)//'.int'
1705       endif
1706       rmsname=prefintin(:lenint)//'.rms'
1707       open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1708              status='unknown')
1709       open (jrms,file=rmsname,status='unknown')
1710       open(iout,file=outname,status='unknown')
1711 ! Get parameter filenames and open the parameter files.
1712       call getenv('BONDPAR',bondname)
1713       open (ibond,file=bondname,status='old')
1714       call getenv('THETPAR',thetname)
1715       open (ithep,file=thetname,status='old')
1716       call getenv('ROTPAR',rotname)
1717       open (irotam,file=rotname,status='old')
1718       call getenv('TORPAR',torname)
1719       open (itorp,file=torname,status='old')
1720       call getenv('TORDPAR',tordname)
1721       open (itordp,file=tordname,status='old')
1722       call getenv('FOURIER',fouriername)
1723       open (ifourier,file=fouriername,status='old')
1724       call getenv('ELEPAR',elename)
1725       open (ielep,file=elename,status='old')
1726       call getenv('SIDEPAR',sidename)
1727       open (isidep,file=sidename,status='old')
1728       call getenv('SIDEP',sidepname)
1729       open (isidep1,file=sidepname,status="old")
1730       call getenv('SCCORPAR',sccorname)
1731       open (isccor,file=sccorname,status="old")
1732       call getenv('THETPAR_NUCL',thetname_nucl)
1733       open (ithep_nucl,file=thetname_nucl,status='old')
1734       call getenv('ROTPAR_NUCL',rotname_nucl)
1735       open (irotam_nucl,file=rotname_nucl,status='old')
1736       call getenv('TORPAR_NUCL',torname_nucl)
1737       open (itorp_nucl,file=torname_nucl,status='old')
1738       call getenv('TORDPAR_NUCL',tordname_nucl)
1739       open (itordp_nucl,file=tordname_nucl,status='old')
1740       call getenv('SIDEPAR_NUCL',sidename_nucl)
1741       open (isidep_nucl,file=sidename_nucl,status='old')
1742       call getenv('SIDEPAR_SCBASE',sidename_scbase)
1743       open (isidep_scbase,file=sidename_scbase,status='old')
1744       call getenv('PEPPAR_PEPBASE',pepname_pepbase)
1745       open (isidep_pepbase,file=pepname_pepbase,status='old')
1746       call getenv('SCPAR_PHOSPH',pepname_scpho)
1747       open (isidep_scpho,file=pepname_scpho,status='old')
1748       call getenv('PEPPAR_PHOSPH',pepname_peppho)
1749       open (isidep_peppho,file=pepname_peppho,status='old')
1750
1751
1752       call getenv('LIPTRANPAR',liptranname)
1753       open (iliptranpar,file=liptranname,status='old')
1754       call getenv('TUBEPAR',tubename)
1755       open (itube,file=tubename,status='old')
1756       call getenv('IONPAR',ionname)
1757       open (iion,file=ionname,status='old')
1758
1759 #ifndef OLDSCP
1760 !
1761 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1762 ! Use -DOLDSCP to use hard-coded constants instead.
1763 !
1764       call getenv('SCPPAR',scpname)
1765       open (iscpp,file=scpname,status='old')
1766 #endif
1767       return
1768       end subroutine openunits
1769 !-----------------------------------------------------------------------------
1770 ! geomout.F
1771 !-----------------------------------------------------------------------------
1772       subroutine pdboutC(etot,rmsd,tytul)
1773
1774       use energy_data, only: ihpb,jhpb,itype,molnum
1775 !      implicit real*8 (a-h,o-z)
1776 !      include 'DIMENSIONS'
1777 !      include 'COMMON.CONTROL'
1778 !      include 'COMMON.CHAIN'
1779 !      include 'COMMON.INTERACT'
1780 !      include 'COMMON.NAMES'
1781 !      include 'COMMON.IOUNITS'
1782 !      include 'COMMON.HEADER'
1783 !      include 'COMMON.SBRIDGE'
1784 !      include 'COMMON.TEMPFAC'
1785       character(len=50) :: tytul
1786       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1787                                                'G','H','I','J'/)
1788       integer :: ica(nres)
1789       real(kind=8) :: etot,rmsd
1790       integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
1791       real(kind=8) :: boxxxx(3)
1792       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1793         ' ENERGY ',etot,' RMS ',rmsd
1794       iatom=0
1795       ichain=1
1796       ires=0
1797       boxxxshift(1)=int(c(1,nnt)/boxxsize)
1798 !      if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
1799       boxxxshift(2)=int(c(2,nnt)/boxzsize)
1800 !      if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
1801       boxxxshift(3)=int(c(3,nnt)/boxzsize)
1802 !      if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
1803
1804       boxxxx(1)=boxxsize
1805       boxxxx(2)=boxysize
1806       boxxxx(3)=boxzsize
1807
1808       do i=nnt,nct
1809         mnum=molnum(i)
1810         iti=itype(i,mnum)
1811         if (iti.eq.ntyp1_molec(mnum)) then
1812           ichain=ichain+1
1813           ires=0
1814           write (ipdb,'(a)') 'TER'
1815         else
1816         ires=ires+1
1817         iatom=iatom+1
1818         ica(i)=iatom
1819         if (mnum.eq.5) then
1820            do j=1,3
1821             if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then
1822            c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j)
1823             else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then
1824            c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j)
1825              else
1826            c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j)
1827             endif
1828            enddo
1829         endif
1830         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
1831            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
1832         if ((iti.ne.10).and.(mnum.ne.5)) then
1833           iatom=iatom+1
1834           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
1835             ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
1836         endif
1837         endif
1838       enddo
1839       write (ipdb,'(a)') 'TER'
1840       do i=nnt,nct-1
1841         mnum=molnum(i)
1842         mnum1=molnum(i+1)
1843         if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
1844         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1845           write (ipdb,30) ica(i),ica(i+1)
1846         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1847           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
1848         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
1849           write (ipdb,30) ica(i),ica(i)+1
1850         endif
1851       enddo
1852       if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
1853         write (ipdb,30) ica(nct),ica(nct)+1
1854       endif
1855       do i=1,nss
1856         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1857       enddo
1858       write (ipdb,'(a6)') 'ENDMDL'
1859   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1860   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1861   30  FORMAT ('CONECT',8I5)
1862       return
1863       end subroutine pdboutC
1864 !-----------------------------------------------------------------------------
1865       subroutine cartout(igr,i,etot,free,rmsd,plik)
1866 !     implicit real*8 (a-h,o-z)
1867 !     include 'DIMENSIONS'
1868 !     include 'sizesclu.dat'
1869 !     include 'COMMON.IOUNITS'
1870 !     include 'COMMON.CHAIN'
1871 !     include 'COMMON.VAR'
1872 !     include 'COMMON.LOCAL'
1873 !     include 'COMMON.INTERACT'
1874 !     include 'COMMON.NAMES'
1875 !     include 'COMMON.GEO'
1876 !     include 'COMMON.CLUSTER'
1877       integer :: igr,i,j,k
1878       real(kind=8) :: etot,free,rmsd
1879       character(len=80) :: plik
1880       open (igeom,file=plik,position='append')
1881       write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
1882       write (igeom,'(i4,$)') &
1883         nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
1884       write (igeom,'(i10)') iscore(i)
1885       write (igeom,'(8f10.5)') &
1886         ((allcart(k,j,i),k=1,3),j=1,nres),&
1887         ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
1888       return
1889       end subroutine cartout
1890 !------------------------------------------------------------------------------
1891 !      subroutine alloc_clust_arrays(n_conf)
1892
1893 !      integer :: n_conf
1894 !COMMON.CLUSTER
1895 !      common /clu/
1896 !      allocate(diss(maxdist)) !(maxdist)
1897 !el      allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
1898 !      allocatable :: enetb !(max_ene,maxstr_proc)
1899 !el      allocate(entfac(maxconf)) !(maxconf)
1900 !      allocatable :: totfree_gr !(maxgr)
1901 !el      allocate(rcutoff(max_cut+1)) !(max_cut+1)
1902 !      common /clu1/
1903 !      allocatable :: licz,iass !(maxgr)
1904 !      allocatable :: nconf !(maxgr,maxingr)
1905 !      allocatable :: iass_tot !(maxgr,max_cut)
1906 !      allocatable :: list_conf !(maxconf)
1907 !      common /alles/
1908 !el      allocatable :: allcart !(3,maxres2,maxstr_proc)
1909 !el      allocate(rmstb(maxconf)) !(maxconf)
1910 !el      allocate(mult(nres)) !(maxres)
1911 !el      allocatable :: nss_all !(maxstr_proc)
1912 !el      allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
1913 !      allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
1914 !COMMON.TEMPFAC
1915 !      common /factemp/
1916 !      allocatable :: tempfac !(2,maxres)
1917 !COMMON.FREE
1918 !      common /free/
1919 !el      allocate(beta_h(maxT)) !(maxT)
1920
1921 !      end subroutine alloc_clust_arrays
1922 !-----------------------------------------------------------------------------
1923 !-----------------------------------------------------------------------------
1924       end module io_clust