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       nss_all=0
 563       ihpb_all=0
 564       jhpb_all=0
 565
 566 #ifdef CHUJ
 567       ICON=1
 568   123 continue
 569       if (from_cart .and. .not. from_bx .and. .not. from_cx) then
 570         if (efree) then
 571         read (intin,*,end=13,err=11) energy(icon),totfree(icon),&
 572           rmstb(icon),&
 573           nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
 574           i=1,nss_all(icon)),iscore(icon)
 575         else
 576         read (intin,*,end=13,err=11) energy(icon),rmstb(icon),&
 577           nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
 578           i=1,nss_all(icon)),iscore(icon)
 579         endif
 580         read (intin,'(8f10.5)',end=13,err=10) &
 581           ((allcart(j,i,icon),j=1,3),i=1,nres),&
 582           ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
 583         print *,icon,energy(icon),nss_all(icon),rmstb(icon)
 584       else
 585         read(intin,'(a80)',end=13,err=12) lineh
 586         read(lineh(:5),*,err=8) ic
 587         if (efree) then
 588         read(lineh(6:),*,err=8) energy(icon)
 589         else
 590         read(lineh(6:),*,err=8) energy(icon)
 591         endif
 592         goto 9
 593     8   ic=1
 594         print *,'error, assuming e=1d10',lineh
 595         energy(icon)=1d10
 596         nss=0
 597     9   continue
 598 !old        read(lineh(18:),*,end=13,err=11) nss_all(icon)
 599         ii = index(lineh(15:)," ")+15
 600         read(lineh(ii:),*,end=13,err=11) nss_all(icon)
 601         IF (NSS_all(icon).LT.9) THEN
 602           read (lineh(20:),*,end=102) &
 603           (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),&
 604           iscore(icon)
 605         ELSE
 606           read (lineh(20:),*,end=102) &
 607                  (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
 608           read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),&
 609             I=9,NSS_all(icon)),iscore(icon)
 610         ENDIF
 611
 612   102   continue
 613
 614         PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
 615         call read_angles(intin,*13)
 616         do i=1,nres
 617           phiall(i,icon)=phi(i)
 618           thetall(i,icon)=theta(i)
 619           alphall(i,icon)=alph(i)
 620           omall(i,icon)=omeg(i)
 621         enddo
 622       endif
 623       ICON=ICON+1
 624       GOTO 123
 625 !
 626 ! CALCULATE DISTANCES
 627 !
 628    10 print *,'something wrong with angles'
 629       goto 13
 630    11 print *,'something wrong with NSS',nss
 631       goto 13
 632    12 print *,'something wrong with header'
 633
 634    13 NCON=ICON-1
 635
 636 #endif
 637       call flush(iout)
 638       jj_old=1
 639       open (icbase,file=bprotfiles,status="unknown",&
 640          form="unformatted",access="direct",recl=lenrec)
 641 ! Read conformations from binary DA files (one per batch) and write them to
 642 ! a binary DA scratchfile.
 643       jj=0
 644       jjj=0
 645 #ifdef MPI
 646       write (liczba,'(bz,i3.3)') me
 647       IF (ME.EQ.MASTER) THEN
 648 ! Only the master reads the database; it'll send it to the other procs
 649 ! through a ring.
 650 #endif
 651         t_acq = tcpu()
 652         icount=0
 653
 654         if (from_bx) then
 655
 656           open (intin,file=intinname,status="old",form="unformatted",&
 657                   access="direct",recl=lenrec_in)
 658
 659         else if (from_cx) then
 660 #if (defined(AIX) && !defined(JUBL))
 661           call xdrfopen_(ixdrf,intinname, "r", iret)
 662 #else
 663           call xdrfopen(ixdrf,intinname, "r", iret)
 664 #endif
 665           prec=10000.0
 666           write (iout,*) "xdrfopen: iret",iret
 667           if (iret.eq.0) then
 668             write (iout,*) "Error: coordinate file ",&
 669              intinname(:ilen(intinname))," does not exist."
 670             call flush(iout)
 671 #ifdef MPI
 672             call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
 673 #endif
 674             stop
 675           endif
 676         else
 677           write (iout,*) "Error: coordinate format not specified"
 678           call flush(iout)
 679 #ifdef MPI
 680           call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
 681 #else
 682           stop
 683 #endif
 684         endif
 685
 686 !#define DEBUG
 687 #ifdef DEBUG
 688         write (iout,*) "Opening file ",intinname(:ilen(intinname))
 689         write (iout,*) "lenrec",lenrec_in
 690         call flush(iout)
 691 #endif
 692 !#undef DEBUG
 693 !        write (iout,*) "maxconf",maxconf
 694         i=0
 695         do while (.true.)
 696            i=i+1
 697 !el           if (i.gt.maxconf) then
 698 !el             write (iout,*) "Error: too many conformations ",&
 699 !el              "(",maxconf,") maximum."
 700 !#ifdef MPI
 701 !el             call MPI_Abort(MPI_COMM_WORLD,errcode,ierror)
 702 !#endif
 703 !el             stop
 704 !el           endif
 705 !          write (iout,*) "i",i
 706 !          call flush(iout)
 707           if (from_bx) then
 708             read(intin,err=101,end=101) &
 709              ((csingle(l,k),l=1,3),k=1,nres),&
 710              ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 711              nss,(ihpb(k),jhpb(k),k=1,nss),&
 712              energy(jj+1),&
 713              entfac(jj+1),rmstb(jj+1),iscor
 714              do j=1,2*nres
 715                do k=1,3
 716                  c(k,j)=csingle(k,j)
 717                enddo
 718              enddo
 719           else
 720 #if (defined(AIX) && !defined(JUBL))
 721             call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
 722             if (iret.eq.0) goto 101
 723             call xdrfint_(ixdrf, nss, iret)
 724             if (iret.eq.0) goto 101
 725             do j=1,nss
 726               call xdrfint_(ixdrf, ihpb(j), iret)
 727               if (iret.eq.0) goto 101
 728               call xdrfint_(ixdrf, jhpb(j), iret)
 729               if (iret.eq.0) goto 101
 730             enddo
 731             call xdrffloat_(ixdrf,reini,iret)
 732             if (iret.eq.0) goto 101
 733             call xdrffloat_(ixdrf,refree,iret)
 734             if (iret.eq.0) goto 101
 735             call xdrffloat_(ixdrf,rmsdev,iret)
 736             if (iret.eq.0) goto 101
 737             call xdrfint_(ixdrf,iscor,iret)
 738             if (iret.eq.0) goto 101
 739 #else
 740 !            write (iout,*) "calling xdrf3dfcoord"
 741             call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
 742 !            write (iout,*) "iret",iret
 743 !            call flush(iout)
 744             if (iret.eq.0) goto 101
 745             call xdrfint(ixdrf, nss, iret)
 746             write (iout,*) "iret",iret
 747 !            write (iout,*) "nss",nss,i,"TUTU"
 748             call flush(iout)
 749             if (iret.eq.0) goto 101
 750             do k=1,nss
 751               call xdrfint(ixdrf, ihpb(k), iret)
 752               if (iret.eq.0) goto 101
 753               call xdrfint(ixdrf, jhpb(k), iret)
 754               if (iret.eq.0) goto 101
 755 !            write(iout,*), ihpb(k),jhpb(k),"TUTU"
 756             enddo
 757             call xdrffloat(ixdrf,reini,iret)
 758             if (iret.eq.0) goto 101
 759 !            write(iout,*) "TUTU", reini
 760             call xdrffloat(ixdrf,refree,iret)
 761 !            write(iout,*) "TUTU", refree
 762             if (iret.eq.0) goto 101
 763             call xdrffloat(ixdrf,rmsdev,iret)
 764             if (iret.eq.0) goto 101
 765             call xdrfint(ixdrf,iscor,iret)
 766             if (iret.eq.0) goto 101
 767 #endif
 768             energy(jj+1)=reini
 769             entfac(jj+1)=refree
 770             rmstb(jj+1)=rmsdev
 771             do k=1,nres
 772               do l=1,3
 773                 c(l,k)=csingle(l,k)
 774               enddo
 775             enddo
 776             do k=nnt,nct
 777               do l=1,3
 778                 c(l,nres+k)=csingle(l,nres+k-nnt+1)
 779               enddo
 780             enddo
 781           endif
 782 #ifdef DEBUG
 783           write (iout,'(5hREAD ,i5,3f15.4,i10)') &
 784            jj+1,energy(jj+1),entfac(jj+1),&
 785            rmstb(jj+1),iscor
 786           write (iout,*) "Conformation",jjj+1,jj+1
 787           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 788           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 789           call flush(iout)
 790 #endif
 791           call add_new_cconf(jjj,jj,jj_old,icount,Next)
 792         enddo
 793   101   continue
 794         write (iout,*) i-1," conformations read from DA file ",&
 795           intinname(:ilen(intinname))
 796         write (iout,*) jj," conformations read so far"
 797         if (from_bx) then
 798           close(intin)
 799         else
 800 #if (defined(AIX) && !defined(JUBL))
 801           call xdrfclose_(ixdrf, iret)
 802 #else
 803           call xdrfclose(ixdrf, iret)
 804 #endif
 805         endif
 806 #ifdef MPI
 807 !#ifdef DEBUG
 808         write (iout,*) "jj_old",jj_old," jj",jj
 809 !#endif
 810         call write_and_send_cconf(icount,jj_old,jj,Next)
 811         call MPI_Send(0,1,MPI_INTEGER,Next,570,&
 812                    MPI_COMM_WORLD,IERROR)
 813         jj_old=jj+1
 814 #else
 815         call write_and_send_cconf(icount,jj_old,jj,Next)
 816 #endif
 817         t_acq = tcpu() - t_acq
 818 #ifdef MPI
 819         write (iout,*) "Processor",me,&
 820           " time for conformation read/send",t_acq
 821       ELSE
 822 ! A worker gets the confs from the master and sends them to its neighbor
 823         t_acq = tcpu()
 824         call receive_and_pass_cconf(icount,jj_old,jj,&
 825           Previous,Next)
 826         t_acq = tcpu() - t_acq
 827       ENDIF
 828 #endif
 829       ncon=jj
 830 !      close(icbase)
 831       close(intin)
 832
 833       write(iout,*)"A total of",ncon," conformations read."
 834
 835       allocate(enetb(1:max_ene,ncon)) !(max_ene,maxstr_proc)
 836 #ifdef MPI
 837 ! Check if everyone has the same number of conformations
 838       call MPI_Allgather(ncon,1,MPI_INTEGER,&
 839         ntot_all(0),1,MPI_INTEGER,MPI_Comm_World,IERROR)
 840       lerr=.false.
 841       do i=0,nprocs-1
 842         if (i.ne.me) then
 843             if (ncon.ne.ntot_all(i)) then
 844               write (iout,*) "Number of conformations at processor",i,&
 845                " differs from that at processor",me,&
 846                ncon,ntot_all(i)
 847               lerr = .true.
 848             endif
 849         endif
 850       enddo
 851       if (lerr) then
 852         write (iout,*)
 853         write (iout,*) "Number of conformations read by processors"
 854         write (iout,*)
 855         do i=0,nprocs-1
 856           write (iout,'(8i10)') i,ntot_all(i)
 857         enddo
 858         write (iout,*) "Calculation terminated."
 859         call flush(iout)
 860         return 1
 861       endif
 862       return
 863 #endif
 864  1111 write(iout,*) "Error opening coordinate file ",&
 865        intinname(:ilen(intinname))
 866       call flush(iout)
 867       return 1
 868       end subroutine read_coords
 869 !------------------------------------------------------------------------------
 870       subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
 871
 872       use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
 873       use energy_data, only: itel,itype,dsc,max_ene,molnum
 874       use control_data, only: symetr
 875       use geometry, only: int_from_cart1
 876 !      implicit none
 877 !      include "DIMENSIONS"
 878 !      include "sizesclu.dat"
 879 !      include "COMMON.CLUSTER"
 880 !      include "COMMON.CONTROL"
 881 !      include "COMMON.CHAIN"
 882 !      include "COMMON.INTERACT"
 883 !      include "COMMON.LOCAL"
 884 !      include "COMMON.IOUNITS"
 885 !      include "COMMON.NAMES"
 886 !      include "COMMON.VAR"
 887 !      include "COMMON.SBRIDGE"
 888 !      include "COMMON.GEO"
 889       integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
 890         nn,nn1,inan,Next,itj,chalen,mnum
 891       real(kind=8) :: etot,energia(0:max_ene)
 892       jjj=jjj+1
 893       chalen=int((nct-nnt+2)/symetr)
 894       call int_from_cart1(.false.)
 895       do j=nnt+1,nct-1
 896         mnum=molnum(j)
 897 !        write (iout,*) "Check atom",j
 898         if (mnum.ne.1) cycle
 899         if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
 900         if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle
 901
 902         if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
 903          if (j.gt.2) then
 904           if (itel(j).ne.0 .and. itel(j-1).ne.0) then
 905           write (iout,*) "Conformation",jjj,jj+1
 906           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
 907        chalen
 908           write (iout,*) "The Cartesian geometry is:"
 909           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 910           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 911           write (iout,*) "The internal geometry is:"
 912           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 913           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 914           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 915           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 916           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 917           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 918           write (iout,*) &
 919             "This conformation WILL NOT be added to the database."
 920           return
 921           endif
 922          endif
 923         endif
 924       enddo
 925       do j=nnt,nct
 926         mnum=molnum(j)
 927         itj=itype(j,mnum)
 928         if (mnum.ne.1) cycle
 929         if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) &
 930             .and.  (vbld(nres+j)-dsc(iabs(itj))) &
 931                                   .gt.2.0d0) then
 932           write (iout,*) "Conformation",jjj,jj+1
 933           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
 934           write (iout,*) "The Cartesian geometry is:"
 935           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 936           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 937           write (iout,*) "The internal geometry is:"
 938           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 939           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 940           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 941           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 942           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 943           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 944           write (iout,*) &
 945             "This conformation WILL NOT be added to the database."
 946           return
 947         endif
 948       enddo
 949       do j=3,nres
 950         if (theta(j).le.0.0d0) then
 951           write (iout,*) &
 952             "Zero theta angle(s) in conformation",jjj,jj+1
 953           write (iout,*) "The Cartesian geometry is:"
 954           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 955           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 956           write (iout,*) "The internal geometry is:"
 957           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 958           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 959           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 960           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 961           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 962           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 963           write (iout,*) &
 964             "This conformation WILL NOT be added to the database."
 965           return
 966         endif
 967         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 968       enddo
 969       jj=jj+1
 970 #ifdef DEBUG
 971       write (iout,*) "Conformation",jjj,jj
 972       write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 973       write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 974       write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 975       write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 976       write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 977       write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 978       write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
 979       write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 980       write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 981       write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 982       write (iout,'(e15.5,16i5)') entfac(icount+1)
 983 !     &        iscore(icount+1,0)
 984 #endif
 985       icount=icount+1
 986       call store_cconf_from_file(jj,icount)
 987       if (icount.eq.maxstr_proc) then
 988 #ifdef DEBUG
 989         write (iout,* ) "jj_old",jj_old," jj",jj
 990 #endif
 991         call write_and_send_cconf(icount,jj_old,jj,Next)
 992         jj_old=jj+1
 993         icount=0
 994       endif
 995       return
 996       end subroutine add_new_cconf
 997 !------------------------------------------------------------------------------
 998       subroutine store_cconf_from_file(jj,icount)
 999
1000       use energy_data, only: ihpb,jhpb
1001 !      implicit none
1002 !      include "DIMENSIONS"
1003 !      include "sizesclu.dat"
1004 !      include "COMMON.CLUSTER"
1005 !      include "COMMON.CHAIN"
1006 !      include "COMMON.SBRIDGE"
1007 !      include "COMMON.INTERACT"
1008 !      include "COMMON.IOUNITS"
1009 !      include "COMMON.VAR"
1010       integer :: i,j,jj,icount
1011 ! Store the conformation that has been read in
1012       do i=1,2*nres
1013         do j=1,3
1014           allcart(j,i,icount)=c(j,i)
1015         enddo
1016       enddo
1017       nss_all(icount)=nss
1018       do i=1,nss
1019         ihpb_all(i,icount)=ihpb(i)
1020         jhpb_all(i,icount)=jhpb(i)
1021       enddo
1022       return
1023       end subroutine store_cconf_from_file
1024 !------------------------------------------------------------------------------
1025       subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1026
1027 !      implicit none
1028 !      include "DIMENSIONS"
1029 !      include "sizesclu.dat"
1030 #ifdef MPI
1031       use MPI_data
1032       include "mpif.h"
1033       integer :: IERROR
1034 !      include "COMMON.MPI"
1035 #endif
1036 !      include "COMMON.CHAIN"
1037 !      include "COMMON.SBRIDGE"
1038 !      include "COMMON.INTERACT"
1039 !      include "COMMON.IOUNITS"
1040 !      include "COMMON.CLUSTER"
1041 !      include "COMMON.VAR"
1042       integer :: icount,jj_old,jj,Next
1043 ! Write the structures to a scratch file
1044 #ifdef MPI
1045 ! Master sends the portion of conformations that have been read in to the neighbor
1046 #ifdef DEBUG
1047       write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1048       call flush(iout)
1049 #endif
1050       call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1051       call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1052           Next,571,MPI_COMM_WORLD,IERROR)
1053       call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1054           Next,572,MPI_COMM_WORLD,IERROR)
1055       call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1056           Next,573,MPI_COMM_WORLD,IERROR)
1057       call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1058           Next,577,MPI_COMM_WORLD,IERROR)
1059       call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1060           Next,579,MPI_COMM_WORLD,IERROR)
1061       call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1062           MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1063 #endif
1064       call dawrite_ccoords(jj_old,jj,icbase)
1065       return
1066       end subroutine write_and_send_cconf
1067 !------------------------------------------------------------------------------
1068 #ifdef MPI
1069       subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1070
1071       use MPI_data
1072 !      implicit none
1073 !      include "DIMENSIONS"
1074 !      include "sizesclu.dat"
1075       include "mpif.h"
1076       integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1077 !      include "COMMON.MPI"
1078 !      include "COMMON.CHAIN"
1079 !      include "COMMON.SBRIDGE"
1080 !      include "COMMON.INTERACT"
1081 !      include "COMMON.IOUNITS"
1082 !      include "COMMON.VAR"
1083 !      include "COMMON.GEO"
1084 !      include "COMMON.CLUSTER"
1085       integer :: i,j,k,icount,jj_old,jj,Previous,Next
1086       icount=1
1087 #ifdef DEBUG
1088       write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1089       call flush(iout)
1090 #endif
1091       do while (icount.gt.0)
1092       call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1093            STATUS,IERROR)
1094       call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1095            IERROR)
1096 #ifdef DEBUG
1097       write (iout,*) "Processor",me," icount",icount
1098 #endif
1099       if (icount.eq.0) return
1100       call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1101           Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1102       call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1103         Next,571,MPI_COMM_WORLD,IERROR)
1104       call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1105           Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1106       call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1107         Next,572,MPI_COMM_WORLD,IERROR)
1108       call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1109           Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1110       call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1111         Next,573,MPI_COMM_WORLD,IERROR)
1112       call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1113         Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1114       call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1115         Next,577,MPI_COMM_WORLD,IERROR)
1116       call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1117         Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1118       call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1119         Next,579,MPI_COMM_WORLD,IERROR)
1120       call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1121         MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1122       call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1123         MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1124       jj=jj_old+icount-1
1125       call dawrite_ccoords(jj_old,jj,icbase)
1126       jj_old=jj+1
1127 #ifdef DEBUG
1128       write (iout,*) "Processor",me," received",icount," conformations"
1129       do i=1,icount
1130         write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1131         write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1132         write (iout,'(e15.5,16i5)') entfac(i)
1133       enddo
1134 #endif
1135       enddo
1136       return
1137       end subroutine receive_and_pass_cconf
1138 #endif
1139 !------------------------------------------------------------------------------
1140       subroutine daread_ccoords(istart_conf,iend_conf)
1141
1142       use energy_data, only: dyn_ss
1143
1144 !      implicit none
1145 !      include "DIMENSIONS"
1146 !      include "sizesclu.dat"
1147 #ifdef MPI
1148       use MPI_data
1149       include "mpif.h"
1150 !      include "COMMON.MPI"
1151 #endif
1152 !      include "COMMON.CHAIN"
1153 !      include "COMMON.CLUSTER"
1154 !      include "COMMON.IOUNITS"
1155 !      include "COMMON.INTERACT"
1156 !      include "COMMON.VAR"
1157 !      include "COMMON.SBRIDGE"
1158 !      include "COMMON.GEO"
1159       integer :: istart_conf,iend_conf
1160       integer :: i,j,ij,ii,iii
1161       integer :: len
1162       character(len=16) :: form,acc
1163       character(len=32) :: nam
1164 !
1165 ! Read conformations off a DA scratchfile.
1166 !
1167 #ifdef DEBUG
1168       write (iout,*) "DAREAD_COORDS"
1169       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1170       inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1171       write (iout,*) "len=",len," form=",form," acc=",acc
1172       write (iout,*) "nam=",nam
1173       call flush(iout)
1174 #endif
1175       do ii=istart_conf,iend_conf
1176         ij = ii - istart_conf + 1
1177         iii=list_conf(ii)
1178 #ifdef DEBUG
1179         write (iout,*) "Reading binary file, record",iii," ii",ii
1180         call flush(iout)
1181 #endif
1182         if (dyn_ss) then
1183         read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
1184          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1185          entfac(ii),rmstb(ii)
1186         else
1187         read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1188          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1189          nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1190          entfac(ii),rmstb(ii)
1191         endif
1192
1193 #ifdef DEBUG
1194         write (iout,*) ii,iii,ij,entfac(ii)
1195         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1196         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1197           i=nnt+nres,nct+nres)
1198         write (iout,'(2e15.5)') entfac(ij)
1199         write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1200           jhpb_all(i,ij),i=1,nss)
1201         call flush(iout)
1202 #endif
1203       enddo
1204       return
1205       end subroutine daread_ccoords
1206 !------------------------------------------------------------------------------
1207       subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1208
1209 !      implicit none
1210 !      include "DIMENSIONS"
1211 !      include "sizesclu.dat"
1212       use energy_data, only: dyn_ss
1213
1214 #ifdef MPI
1215       use MPI_data
1216       include "mpif.h"
1217 !      include "COMMON.MPI"
1218 #endif
1219 !      include "COMMON.CHAIN"
1220 !      include "COMMON.INTERACT"
1221 !      include "COMMON.IOUNITS"
1222 !      include "COMMON.VAR"
1223 !      include "COMMON.SBRIDGE"
1224 !      include "COMMON.GEO"
1225 !      include "COMMON.CLUSTER"
1226       integer :: istart_conf,iend_conf
1227       integer :: i,j,ii,ij,iii,unit_out
1228       integer :: len
1229       character(len=16) :: form,acc
1230       character(len=32) :: nam
1231 !
1232 ! Write conformations to a DA scratchfile.
1233 !
1234 #ifdef DEBUG
1235       write (iout,*) "DAWRITE_COORDS"
1236       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1237       write (iout,*) "lenrec",lenrec
1238       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1239       write (iout,*) "len=",len," form=",form," acc=",acc
1240       write (iout,*) "nam=",nam
1241       call flush(iout)
1242 #endif
1243       do ii=istart_conf,iend_conf
1244         iii=list_conf(ii)
1245         ij = ii - istart_conf + 1
1246 #ifdef DEBUG
1247         write (iout,*) "Writing binary file, record",iii," ii",ii
1248         call flush(iout)
1249 #endif
1250        if (dyn_ss) then
1251         write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1252          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1253          entfac(ii),rmstb(ii)
1254         else
1255         write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
1256          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1257          nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1258          entfac(ii),rmstb(ii)
1259        endif
1260
1261 #ifdef DEBUG
1262         write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1263         write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1264          nct+nres)
1265         write (iout,'(2e15.5)') entfac(ij)
1266         write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1267          nss_all(ij))
1268         call flush(iout)
1269 #endif
1270       enddo
1271       return
1272       end subroutine dawrite_ccoords
1273 !-----------------------------------------------------------------------------
1274 ! readrtns.F
1275 !-----------------------------------------------------------------------------
1276       subroutine read_control
1277 !
1278 ! Read molecular data
1279 !
1280       use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
1281       use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1282                  iscode,symetr,punch_dist,print_dist,nstart,nend,&
1283                  caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
1284                  r_cut_ele,nclust,tor_mode,scelemode
1285 !      implicit none
1286 !      include 'DIMENSIONS'
1287 !      include 'sizesclu.dat'
1288 !      include 'COMMON.IOUNITS'
1289 !      include 'COMMON.TIME1'
1290 !      include 'COMMON.SBRIDGE'
1291 !      include 'COMMON.CONTROL'
1292 !      include 'COMMON.CLUSTER'
1293 !      include 'COMMON.CHAIN'
1294 !      include 'COMMON.HEADER'
1295 !      include 'COMMON.FFIELD'
1296 !      include 'COMMON.FREE'
1297 !      include 'COMMON.INTERACT'
1298       character(len=320) :: controlcard !,ucase
1299 !#ifdef MPL
1300 !      include 'COMMON.INFO'
1301 !#endif
1302       integer :: i
1303
1304       read (INP,'(a80)') titel
1305       call card_concat(controlcard,.true.)
1306
1307       call readi(controlcard,'NRES',nres_molec(1),0)
1308       call readi(controlcard,'NRES_NUCL',nres_molec(2),0)
1309       call readi(controlcard,'NRES_CAT',nres_molec(5),0)
1310       nres=0
1311       do i=1,5
1312        nres=nres_molec(i)+nres
1313       enddo
1314       print *,"TU",nres_molec(:)
1315
1316 !      call alloc_clust_arrays
1317       allocate(rcutoff(max_cut+1)) !(max_cut+1)
1318       allocate(beta_h(maxT)) !(maxT)
1319       allocate(mult(nres)) !(maxres)
1320
1321
1322       call readi(controlcard,'RESCALE',rescale_mode,2)
1323       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1324       write (iout,*) "DISTCHAINMAX",distchainmax
1325       call readi(controlcard,'PDBOUT',outpdb,0)
1326       call readi(controlcard,'MOL2OUT',outmol2,0)
1327       refstr=(index(controlcard,'REFSTR').gt.0)
1328       write (iout,*) "REFSTR",refstr
1329       pdbref=(index(controlcard,'PDBREF').gt.0)
1330       iscode=index(controlcard,'ONE_LETTER')
1331       tree=(index(controlcard,'MAKE_TREE').gt.0)
1332       min_var=(index(controlcard,'MINVAR').gt.0)
1333       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1334       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1335       call readi(controlcard,'NCUT',ncut,1)
1336       call readi(controlcard,'SYM',symetr,1)
1337       write (iout,*) 'sym', symetr
1338       call readi(controlcard,'NSTART',nstart,0)
1339       call reada(controlcard,'BOXX',boxxsize,100.0d0)
1340       call reada(controlcard,'BOXY',boxysize,100.0d0)
1341       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1342       call readi(controlcard,'NCLUST',nclust,5)
1343 !      ions=index(controlcard,"IONS").gt.0
1344       call readi(controlcard,"SCELEMODE",scelemode,0)
1345       print *,"SCELE",scelemode
1346       call readi(controlcard,'TORMODE',tor_mode,0)
1347 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1348         write(iout,*) "torsional and valence angle mode",tor_mode
1349       call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
1350       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1351       write(iout,*) "R_CUT_ELE=",r_cut_ele
1352       call readi(controlcard,'NEND',nend,0)
1353       call reada(controlcard,'ECUT',ecut,10.0d0)
1354       call reada(controlcard,'PROB',prob_limit,0.99d0)
1355       write (iout,*) "Probability limit",prob_limit
1356       lgrp=(index(controlcard,'LGRP').gt.0)
1357       caonly=(index(controlcard,'CA_ONLY').gt.0)
1358       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1359       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1360       call readi(controlcard,'IOPT',iopt,2)
1361       lside = index(controlcard,"SIDE").gt.0
1362       efree = index(controlcard,"EFREE").gt.0
1363       call readi(controlcard,'NTEMP',nT,1)
1364       write (iout,*) "nT",nT
1365 !el      call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1366       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1367       write (iout,*) "nT",nT
1368       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1369       do i=1,nT
1370         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1371       enddo
1372       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1373       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1374       lprint_int=index(controlcard,"PRINT_INT") .gt.0
1375       if (min_var) iopt=1
1376       return
1377       end subroutine read_control
1378 !-----------------------------------------------------------------------------
1379       subroutine molread
1380 !
1381 ! Read molecular data.
1382 !
1383       use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
1384       use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1385 !                 wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1386 !                 wturn3,wturn4,wturn6,wvdwpp,weights
1387       use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1388                  indpdb
1389       use geometry, only: chainbuild,alloc_geo_arrays
1390       use energy, only: alloc_ener_arrays
1391       use io_base, only: rescode
1392       use control, only: setup_var,init_int_table
1393       use conform_compar, only: contact
1394 !      implicit none
1395 !      include 'DIMENSIONS'
1396 !      include 'COMMON.IOUNITS'
1397 !      include 'COMMON.GEO'
1398 !      include 'COMMON.VAR'
1399 !      include 'COMMON.INTERACT'
1400 !      include 'COMMON.LOCAL'
1401 !      include 'COMMON.NAMES'
1402 !      include 'COMMON.CHAIN'
1403 !      include 'COMMON.FFIELD'
1404 !      include 'COMMON.SBRIDGE'
1405 !      include 'COMMON.HEADER'
1406 !      include 'COMMON.CONTROL'
1407 !      include 'COMMON.CONTACTS'
1408 !      include 'COMMON.TIME1'
1409 !#ifdef MPL
1410 !      include 'COMMON.INFO'
1411 !#endif
1412       character(len=4) ::  sequence(nres) !(maxres)
1413       character(len=800) :: weightcard
1414 !      integer rescode
1415       real(kind=8) :: x(6*nres) !(maxvar)
1416       real(kind=8) :: ssscale
1417       integer  :: itype_pdb(nres) !(maxres)
1418 !      logical seq_comp
1419       integer :: i,j,kkk,mnum
1420 !
1421 ! Body
1422 !
1423 !el      allocate(weights(n_ene))
1424       allocate(weights(max_ene))
1425       call alloc_geo_arrays
1426       call alloc_ener_arrays
1427 !-----------------------------
1428       allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1429       allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1430       allocate(itype(nres+2,5)) !(maxres)
1431       allocate(molnum(nres+1))
1432       allocate(itel(nres+2))
1433
1434       do i=1,2*nres+2
1435         do j=1,3
1436           c(j,i)=0
1437           dc(j,i)=0
1438         enddo
1439       enddo
1440       do j=1,5
1441       do i=1,nres+2
1442         itype(i,j)=0
1443         itel(i)=0
1444       enddo
1445       enddo
1446 !--------------------------
1447 ! Read weights of the subsequent energy terms.
1448       call card_concat(weightcard,.true.)
1449       call reada(weightcard,'WSC',wsc,1.0d0)
1450       call reada(weightcard,'WLONG',wsc,wsc)
1451       call reada(weightcard,'WSCP',wscp,1.0d0)
1452       call reada(weightcard,'WELEC',welec,1.0D0)
1453       call reada(weightcard,'WVDWPP',wvdwpp,welec)
1454       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1455       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1456       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1457       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1458       call reada(weightcard,'WTURN3',wturn3,1.0D0)
1459       call reada(weightcard,'WTURN4',wturn4,1.0D0)
1460       call reada(weightcard,'WTURN6',wturn6,1.0D0)
1461       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1462       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1463       call reada(weightcard,'WBOND',wbond,1.0D0)
1464       call reada(weightcard,'WTOR',wtor,1.0D0)
1465       call reada(weightcard,'WTORD',wtor_d,1.0D0)
1466       call reada(weightcard,'WANG',wang,1.0D0)
1467       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1468       call reada(weightcard,'SCAL14',scal14,0.4D0)
1469       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1470       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1471       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1472       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
1473       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
1474       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
1475       if (index(weightcard,'SOFT').gt.0) ipot=6
1476       call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
1477       call reada(weightcard,'WELPP',welpp,0.0d0)
1478       call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
1479       call reada(weightcard,'WELPSB',welpsb,0.0D0)
1480       call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
1481       call reada(weightcard,'WELSB',welsb,0.0D0)
1482       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
1483       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
1484       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
1485       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
1486 !      print *,"KUR..",wtor_nucl
1487       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
1488       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
1489       call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
1490 ! 12/1/95 Added weight for the multi-body term WCORR
1491       call reada(weightcard,'WCORRH',wcorr,1.0D0)
1492       call reada(weightcard,'WSCBASE',wscbase,0.0D0)
1493       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
1494       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
1495       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
1496       call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
1497
1498       call reada(weightcard,"D0CM",d0cm,3.78d0)
1499       call reada(weightcard,"AKCM",akcm,15.1d0)
1500       call reada(weightcard,"AKTH",akth,11.0d0)
1501       call reada(weightcard,"AKCT",akct,12.0d0)
1502       call reada(weightcard,"V1SS",v1ss,-1.08d0)
1503       call reada(weightcard,"V2SS",v2ss,7.61d0)
1504       call reada(weightcard,"V3SS",v3ss,13.7d0)
1505       call reada(weightcard,"EBR",ebr,-5.50D0)
1506       call reada(weightcard,"ATRISS",atriss,0.301D0)
1507       call reada(weightcard,"BTRISS",btriss,0.021D0)
1508       call reada(weightcard,"CTRISS",ctriss,1.001D0)
1509       call reada(weightcard,"DTRISS",dtriss,1.001D0)
1510       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
1511       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
1512
1513       call reada(weightcard,"HT",Ht,0.0D0)
1514       if (dyn_ss) then
1515        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
1516         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
1517         akcm=akcm/wsc*ssscale
1518         akth=akth/wsc*ssscale
1519         akct=akct/wsc*ssscale
1520         v1ss=v1ss/wsc*ssscale
1521         v2ss=v2ss/wsc*ssscale
1522         v3ss=v3ss/wsc*ssscale
1523       else
1524         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
1525       endif
1526
1527       if (wcorr4.gt.0.0d0) wcorr=wcorr4
1528       weights(1)=wsc
1529       weights(2)=wscp
1530       weights(3)=welec
1531       weights(4)=wcorr
1532       weights(5)=wcorr5
1533       weights(6)=wcorr6
1534       weights(7)=wel_loc
1535       weights(8)=wturn3
1536       weights(9)=wturn4
1537       weights(10)=wturn6
1538       weights(11)=wang
1539       weights(12)=wscloc
1540       weights(13)=wtor
1541       weights(14)=wtor_d
1542       weights(15)=wstrain
1543       weights(16)=wvdwpp
1544       weights(17)=wbond
1545       weights(18)=scal14
1546 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
1547       weights(21)=wsccor
1548           weights(26)=wvdwpp_nucl
1549           weights(27)=welpp
1550           weights(28)=wvdwpsb
1551           weights(29)=welpsb
1552           weights(30)=wvdwsb
1553           weights(31)=welsb
1554           weights(32)=wbond_nucl
1555           weights(33)=wang_nucl
1556           weights(34)=wsbloc
1557           weights(35)=wtor_nucl
1558           weights(36)=wtor_d_nucl
1559           weights(37)=wcorr_nucl
1560           weights(38)=wcorr3_nucl
1561           weights(41)=wcatcat
1562           weights(42)=wcatprot
1563           weights(46)=wscbase
1564           weights(47)=wscpho
1565           weights(48)=wpeppho
1566           weights(49)=wpeppho
1567           weights(50)=wcatnucl
1568
1569
1570       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1571         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1572         wturn4,wturn6,wsccor
1573    10 format (/'Energy-term weights (unscaled):'// &
1574        'WSCC=   ',f10.6,' (SC-SC)'/ &
1575        'WSCP=   ',f10.6,' (SC-p)'/ &
1576        'WELEC=  ',f10.6,' (p-p electr)'/ &
1577        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1578        'WBOND=  ',f10.6,' (stretching)'/ &
1579        'WANG=   ',f10.6,' (bending)'/ &
1580        'WSCLOC= ',f10.6,' (SC local)'/ &
1581        'WTOR=   ',f10.6,' (torsional)'/ &
1582        'WTORD=  ',f10.6,' (double torsional)'/ &
1583        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1584        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1585        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1586        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1587        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1588        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1589        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1590        'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1591        'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1592
1593       if (wcorr4.gt.0.0d0) then
1594         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1595          'between contact pairs of peptide groups'
1596         write (iout,'(2(a,f5.3/))') &
1597         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1598         'Range of quenching the correlation terms:',2*delt_corr
1599       else if (wcorr.gt.0.0d0) then
1600         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1601          'between contact pairs of peptide groups'
1602       endif
1603       write (iout,'(a,f8.3)') &
1604         'Scaling factor of 1,4 SC-p interactions:',scal14
1605       write (iout,'(a,f8.3)') &
1606         'General scaling factor of SC-p interactions:',scalscp
1607       r0_corr=cutoff_corr-delt_corr
1608       do i=1,20
1609         aad(i,1)=scalscp*aad(i,1)
1610         aad(i,2)=scalscp*aad(i,2)
1611         bad(i,1)=scalscp*bad(i,1)
1612         bad(i,2)=scalscp*bad(i,2)
1613       enddo
1614
1615       call flush(iout)
1616       print *,'indpdb=',indpdb,' pdbref=',pdbref
1617
1618 ! Read sequence if not taken from the pdb file.
1619       if (iscode.gt.0) then
1620         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1621       else
1622         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1623       endif
1624       print *,nres_molec(:),nres
1625 ! Convert sequence to numeric code
1626       do i=1,nres_molec(1)
1627         mnum=1
1628         molnum(i)=1
1629         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1630       enddo
1631       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
1632         mnum=2
1633         molnum(i)=2
1634         write (iout,*),i,sequence(i)
1635         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1636
1637       enddo
1638       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
1639         mnum=5
1640         molnum(i)=5
1641         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1642       enddo
1643       print *,nres
1644       print '(20i4)',(itype(i,molnum(i)),i=1,nres)
1645
1646       do i=1,nres-1
1647 #ifdef PROCOR
1648         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
1649             .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
1650 #else
1651         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
1652 #endif
1653           itel(i)=0
1654 #ifdef PROCOR
1655         else if (iabs(itype(i+1,1)).ne.20) then
1656 #else
1657         else if (iabs(itype(i,1)).ne.20) then
1658 #endif
1659           itel(i)=1
1660         else
1661           itel(i)=2
1662         endif
1663       enddo
1664       write (iout,*) "ITEL"
1665       do i=1,nres-1
1666         write (iout,*) i,itype(i,molnum(i)),itel(i)
1667       enddo
1668
1669       print *,'Call Read_Bridge.'
1670       call read_bridge
1671       nnt=1
1672       nct=nres
1673       print *,'NNT=',NNT,' NCT=',NCT
1674       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1675       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1676       if (nstart.lt.nnt) nstart=nnt
1677       if (nend.gt.nct .or. nend.eq.0) nend=nct
1678       write (iout,*) "nstart",nstart," nend",nend
1679       nres0=nres
1680 !      if (pdbref) then
1681 !        read(inp,'(a)') pdbfile
1682 !        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1683 !        open(ipdbin,file=pdbfile,status='old',err=33)
1684 !        goto 34
1685 !  33    write (iout,'(a)') 'Error opening PDB file.'
1686 !        stop
1687 !  34    continue
1688 !        print *,'Begin reading pdb data'
1689 !        call readpdb
1690 !        print *,'Finished reading pdb data'
1691 !        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1692 !        do i=1,nres
1693 !          itype_pdb(i)=itype(i)
1694 !        enddo
1695 !        close (ipdbin)
1696 !        write (iout,'(a,i3)') 'nsup=',nsup
1697 !        nstart_seq=nnt
1698 !        if (nsup.le.(nct-nnt+1)) then
1699 !          do i=0,nct-nnt+1-nsup
1700 !            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1701 !              nstart_seq=nnt+i
1702 !              goto 111
1703 !            endif
1704 !          enddo
1705 !          write (iout,'(a)')
1706 !     &            'Error - sequences to be superposed do not match.'
1707 !          stop
1708 !        else
1709 !          do i=0,nsup-(nct-nnt+1)
1710 !            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1711 !     &      then
1712 !              nstart_sup=nstart_sup+i
1713 !              nsup=nct-nnt+1
1714 !              goto 111
1715 !            endif
1716 !          enddo
1717 !          write (iout,'(a)')
1718 !     &            'Error - sequences to be superposed do not match.'
1719 !        endif
1720 !  111   continue
1721 !        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1722 !     &                 ' nstart_seq=',nstart_seq
1723 !      endif
1724 write(iout,*)"przed ini_int_tab"
1725       call init_int_table
1726 write(iout,*)"po ini_int_tab"
1727 write(iout,*)"przed setup var"
1728       call setup_var
1729 write(iout,*)"po setup var"
1730       if (ns.gt.0.and.dyn_ss) then
1731           do i=nss+1,nhpb
1732             ihpb(i-nss)=ihpb(i)
1733             jhpb(i-nss)=jhpb(i)
1734             forcon(i-nss)=forcon(i)
1735             dhpb(i-nss)=dhpb(i)
1736           enddo
1737           nhpb=nhpb-nss
1738           nss=0
1739           link_start=1
1740           link_end=nhpb
1741 !          call hpb_partition
1742           do i=1,ns
1743             dyn_ss_mask(iss(i))=.true.
1744           enddo
1745       endif
1746
1747       write (iout,*) "molread: REFSTR",refstr
1748       if (refstr) then
1749         if (.not.pdbref) then
1750           call read_angles(inp,*38)
1751           goto 39
1752    38     write (iout,'(a)') 'Error reading reference structure.'
1753 #ifdef MPL
1754           call mp_stopall(Error_Msg)
1755 #else
1756           stop 'Error reading reference structure'
1757 #endif
1758    39     call chainbuild
1759           nstart_sup=nnt
1760           nstart_seq=nnt
1761           nsup=nct-nnt+1
1762           kkk=1
1763           do i=1,2*nres
1764             do j=1,3
1765               cref(j,i,kkk)=c(j,i)
1766             enddo
1767           enddo
1768         endif
1769         call contact(.true.,ncont_ref,icont_ref)
1770       endif
1771       return
1772       end subroutine molread
1773 !-----------------------------------------------------------------------------
1774       subroutine openunits
1775 !      implicit none
1776 !      include 'DIMENSIONS'
1777       use control_data, only: from_cx,from_bx,from_cart
1778 #ifdef MPI
1779       use MPI_data
1780       include "mpif.h"
1781       character(len=3) :: liczba
1782 !      include "COMMON.MPI"
1783 #endif
1784 !      include 'COMMON.IOUNITS'
1785 !      include 'COMMON.CONTROL'
1786       integer :: lenpre,lenpot !,ilen
1787 !      external ilen
1788       character(len=16) :: cformat,cprint
1789 !      character(len=16) ucase
1790       integer :: lenint,lenout
1791       call getenv('INPUT',prefix)
1792       call getenv('OUTPUT',prefout)
1793       call getenv('INTIN',prefintin)
1794       call getenv('COORD',cformat)
1795       call getenv('PRINTCOOR',cprint)
1796       call getenv('SCRATCHDIR',scratchdir)
1797       from_bx=.true.
1798       from_cx=.false.
1799       if (index(ucase(cformat),'CX').gt.0) then
1800         from_cx=.true.
1801         from_bx=.false.
1802       endif
1803       from_cart=.true.
1804       lenpre=ilen(prefix)
1805       lenout=ilen(prefout)
1806       lenint=ilen(prefintin)
1807 ! Get the names and open the input files
1808       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1809 #ifdef MPI
1810       write (liczba,'(bz,i3.3)') me
1811       outname=prefout(:lenout)//'_clust.out_'//liczba
1812 #else
1813       outname=prefout(:lenout)//'_clust.out'
1814 #endif
1815       if (from_bx) then
1816         intinname=prefintin(:lenint)//'.bx'
1817       else if (from_cx) then
1818         intinname=prefintin(:lenint)//'.cx'
1819       else
1820         intinname=prefintin(:lenint)//'.int'
1821       endif
1822       rmsname=prefintin(:lenint)//'.rms'
1823       open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1824              status='unknown')
1825       open (jrms,file=rmsname,status='unknown')
1826       open(iout,file=outname,status='unknown')
1827 ! Get parameter filenames and open the parameter files.
1828       call getenv('BONDPAR',bondname)
1829       open (ibond,file=bondname,status='old')
1830       call getenv('THETPAR',thetname)
1831       open (ithep,file=thetname,status='old')
1832       call getenv('ROTPAR',rotname)
1833       open (irotam,file=rotname,status='old')
1834       call getenv('TORPAR',torname)
1835       open (itorp,file=torname,status='old')
1836       call getenv('TORDPAR',tordname)
1837       open (itordp,file=tordname,status='old')
1838       call getenv('FOURIER',fouriername)
1839       open (ifourier,file=fouriername,status='old')
1840       call getenv('ELEPAR',elename)
1841       open (ielep,file=elename,status='old')
1842       call getenv('SIDEPAR',sidename)
1843       open (isidep,file=sidename,status='old')
1844       call getenv('SIDEP',sidepname)
1845       open (isidep1,file=sidepname,status="old")
1846       call getenv('SCCORPAR',sccorname)
1847       open (isccor,file=sccorname,status="old")
1848       call getenv('BONDPAR_NUCL',bondname_nucl)
1849       open (ibond_nucl,file=bondname_nucl,status='old')
1850       call getenv('THETPAR_NUCL',thetname_nucl)
1851       open (ithep_nucl,file=thetname_nucl,status='old')
1852       call getenv('ROTPAR_NUCL',rotname_nucl)
1853       open (irotam_nucl,file=rotname_nucl,status='old')
1854       call getenv('TORPAR_NUCL',torname_nucl)
1855       open (itorp_nucl,file=torname_nucl,status='old')
1856       call getenv('TORDPAR_NUCL',tordname_nucl)
1857       open (itordp_nucl,file=tordname_nucl,status='old')
1858       call getenv('SIDEPAR_NUCL',sidename_nucl)
1859       open (isidep_nucl,file=sidename_nucl,status='old')
1860       call getenv('SIDEPAR_SCBASE',sidename_scbase)
1861       open (isidep_scbase,file=sidename_scbase,status='old')
1862       call getenv('PEPPAR_PEPBASE',pepname_pepbase)
1863       open (isidep_pepbase,file=pepname_pepbase,status='old')
1864       call getenv('SCPAR_PHOSPH',pepname_scpho)
1865       open (isidep_scpho,file=pepname_scpho,status='old')
1866       call getenv('PEPPAR_PHOSPH',pepname_peppho)
1867       open (isidep_peppho,file=pepname_peppho,status='old')
1868
1869
1870       call getenv('LIPTRANPAR',liptranname)
1871       open (iliptranpar,file=liptranname,status='old')
1872       call getenv('TUBEPAR',tubename)
1873       open (itube,file=tubename,status='old')
1874       call getenv('IONPAR',ionname)
1875       open (iion,file=ionname,status='old')
1876       call getenv('IONPAR_NUCL',ionnuclname)
1877       open (iionnucl,file=ionnuclname,status='old')
1878
1879 #ifndef OLDSCP
1880 !
1881 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1882 ! Use -DOLDSCP to use hard-coded constants instead.
1883 !
1884       call getenv('SCPPAR',scpname)
1885       open (iscpp,file=scpname,status='old')
1886       call getenv('SCPPAR_NUCL',scpname_nucl)
1887       open (iscpp_nucl,file=scpname_nucl,status='old')
1888
1889 #endif
1890       return
1891       end subroutine openunits
1892 !-----------------------------------------------------------------------------
1893 ! geomout.F
1894 !-----------------------------------------------------------------------------
1895       subroutine pdboutC(etot,rmsd,tytul)
1896
1897       use energy_data, only: ihpb,jhpb,itype,molnum
1898       use energy, only:boxshift,to_box
1899 !      implicit real*8 (a-h,o-z)
1900 !      include 'DIMENSIONS'
1901 !      include 'COMMON.CONTROL'
1902 !      include 'COMMON.CHAIN'
1903 !      include 'COMMON.INTERACT'
1904 !      include 'COMMON.NAMES'
1905 !      include 'COMMON.IOUNITS'
1906 !      include 'COMMON.HEADER'
1907 !      include 'COMMON.SBRIDGE'
1908 !      include 'COMMON.TEMPFAC'
1909       character(len=50) :: tytul
1910       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1911                                                'G','H','I','J'/)
1912       integer :: ica(nres),k,iti1
1913       real(kind=8) :: etot,rmsd
1914       integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
1915       real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
1916       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1917         ' ENERGY ',etot,' RMS ',rmsd
1918       iatom=0
1919       ichain=1
1920       ires=0
1921 !      boxxxshift(1)=int(c(1,nnt)/boxxsize)
1922 !      if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
1923 !      boxxxshift(2)=int(c(2,nnt)/boxzsize)
1924 !      if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
1925 !      boxxxshift(3)=int(c(3,nnt)/boxzsize)
1926 !      if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
1927 !      do i=1,3
1928 !        if (boxxxshift(i).gt.2) boxxxshift(i)=2
1929 !        if (boxxxshift(i).lt.-2) boxxxshift(i)=-2
1930 !
1931 !      enddo
1932       do j=1,3
1933       cbeg(j)=c(j,nnt)
1934       enddo
1935       call to_box(cbeg(1),cbeg(2),cbeg(3))
1936       boxxxx(1)=boxxsize
1937       boxxxx(2)=boxysize
1938       boxxxx(3)=boxzsize
1939
1940       do i=nnt,nct
1941         mnum=molnum(i)
1942         iti=itype(i,mnum)
1943         if (i.ne.nct) then
1944         iti1=itype(i+1,mnum)
1945         mnum1=molnum(i+1)
1946         if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle
1947         endif
1948         if (iti.eq.ntyp1_molec(mnum)) then
1949           ichain=ichain+1
1950           ires=0
1951           write (ipdb,'(a)') 'TER'
1952         else
1953         ires=ires+1
1954         if ((ires.eq.2).and.(mnum.ne.5)) then
1955          do j=1,3
1956          cbeg(j)=c(j,i-1)
1957          enddo
1958         endif
1959         iatom=iatom+1
1960         ica(i)=iatom
1961         call to_box(c(1,i),c(2,i),c(3,i))
1962
1963         DO k=1,3
1964         Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k))
1965         c(k,i)=cbeg(k)+Rdist(k)
1966         ENDDO
1967       if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then
1968         call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres))
1969         DO k=1,3
1970         Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k))
1971         c(k,i+nres)=cbeg(k)+Rdist(k)
1972         ENDDO
1973        endif
1974         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
1975            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
1976         if ((iti.ne.10).and.(mnum.ne.5)) then
1977           iatom=iatom+1
1978           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
1979             ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
1980         endif
1981         endif
1982       enddo
1983       write (ipdb,'(a)') 'TER'
1984       do i=nnt,nct-1
1985         mnum=molnum(i)
1986         mnum1=molnum(i+1)
1987         if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
1988         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1989           write (ipdb,30) ica(i),ica(i+1)
1990         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1991           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
1992         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
1993           write (ipdb,30) ica(i),ica(i)+1
1994         endif
1995       enddo
1996       if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
1997         write (ipdb,30) ica(nct),ica(nct)+1
1998       endif
1999       do i=1,nss
2000         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2001       enddo
2002       write (ipdb,'(a6)') 'ENDMDL'
2003   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2004   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2005   30  FORMAT ('CONECT',8I5)
2006       return
2007       end subroutine pdboutC
2008 !-----------------------------------------------------------------------------
2009       subroutine cartout(igr,i,etot,free,rmsd,plik)
2010 !     implicit real*8 (a-h,o-z)
2011 !     include 'DIMENSIONS'
2012 !     include 'sizesclu.dat'
2013 !     include 'COMMON.IOUNITS'
2014 !     include 'COMMON.CHAIN'
2015 !     include 'COMMON.VAR'
2016 !     include 'COMMON.LOCAL'
2017 !     include 'COMMON.INTERACT'
2018 !     include 'COMMON.NAMES'
2019 !     include 'COMMON.GEO'
2020 !     include 'COMMON.CLUSTER'
2021       integer :: igr,i,j,k
2022       real(kind=8) :: etot,free,rmsd
2023       character(len=80) :: plik
2024       open (igeom,file=plik,position='append')
2025       write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
2026       write (igeom,'(i4,$)') &
2027         nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
2028       write (igeom,'(i10)') iscore(i)
2029       write (igeom,'(8f10.5)') &
2030         ((allcart(k,j,i),k=1,3),j=1,nres),&
2031         ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
2032       return
2033       end subroutine cartout
2034 !------------------------------------------------------------------------------
2035 !      subroutine alloc_clust_arrays(n_conf)
2036
2037 !      integer :: n_conf
2038 !COMMON.CLUSTER
2039 !      common /clu/
2040 !      allocate(diss(maxdist)) !(maxdist)
2041 !el      allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
2042 !      allocatable :: enetb !(max_ene,maxstr_proc)
2043 !el      allocate(entfac(maxconf)) !(maxconf)
2044 !      allocatable :: totfree_gr !(maxgr)
2045 !el      allocate(rcutoff(max_cut+1)) !(max_cut+1)
2046 !      common /clu1/
2047 !      allocatable :: licz,iass !(maxgr)
2048 !      allocatable :: nconf !(maxgr,maxingr)
2049 !      allocatable :: iass_tot !(maxgr,max_cut)
2050 !      allocatable :: list_conf !(maxconf)
2051 !      common /alles/
2052 !el      allocatable :: allcart !(3,maxres2,maxstr_proc)
2053 !el      allocate(rmstb(maxconf)) !(maxconf)
2054 !el      allocate(mult(nres)) !(maxres)
2055 !el      allocatable :: nss_all !(maxstr_proc)
2056 !el      allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
2057 !      allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
2058 !COMMON.TEMPFAC
2059 !      common /factemp/
2060 !      allocatable :: tempfac !(2,maxres)
2061 !COMMON.FREE
2062 !      common /free/
2063 !el      allocate(beta_h(maxT)) !(maxT)
2064
2065 !      end subroutine alloc_clust_arrays
2066 !-----------------------------------------------------------------------------
2067 !-----------------------------------------------------------------------------
2068       end module io_clust