source/cluster/io_clust.f90

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