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,itype(j,mnum),vbld(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,itype(j-1,molnum(j-1)),itype(j,mnum)
 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,r_cut_mart,r_cut_ang,&
1285                  rlamb_mart
1286       use geometry_data, only:bordliptop,bordlipbot,&
1287           bufliptop,buflipbot,lipthick,lipbufthick
1288 !      implicit none
1289 !      include 'DIMENSIONS'
1290 !      include 'sizesclu.dat'
1291 !      include 'COMMON.IOUNITS'
1292 !      include 'COMMON.TIME1'
1293 !      include 'COMMON.SBRIDGE'
1294 !      include 'COMMON.CONTROL'
1295 !      include 'COMMON.CLUSTER'
1296 !      include 'COMMON.CHAIN'
1297 !      include 'COMMON.HEADER'
1298 !      include 'COMMON.FFIELD'
1299 !      include 'COMMON.FREE'
1300 !      include 'COMMON.INTERACT'
1301       character(len=320) :: controlcard !,ucase
1302 !#ifdef MPL
1303 !      include 'COMMON.INFO'
1304 !#endif
1305       integer :: i
1306
1307       read (INP,'(a80)') titel
1308       call card_concat(controlcard,.true.)
1309
1310       call readi(controlcard,'NRES',nres_molec(1),0)
1311       call readi(controlcard,'NRES_NUCL',nres_molec(2),0)
1312       call readi(controlcard,'NRES_CAT',nres_molec(5),0)
1313       nres=0
1314       do i=1,5
1315        nres=nres_molec(i)+nres
1316       enddo
1317       print *,"TU",nres_molec(:)
1318
1319 !      call alloc_clust_arrays
1320       allocate(rcutoff(max_cut+1)) !(max_cut+1)
1321       allocate(beta_h(maxT)) !(maxT)
1322       allocate(mult(nres)) !(maxres)
1323
1324
1325       call readi(controlcard,'RESCALE',rescale_mode,2)
1326       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1327       write (iout,*) "DISTCHAINMAX",distchainmax
1328       call readi(controlcard,'PDBOUT',outpdb,0)
1329       call readi(controlcard,'MOL2OUT',outmol2,0)
1330       refstr=(index(controlcard,'REFSTR').gt.0)
1331       write (iout,*) "REFSTR",refstr
1332       pdbref=(index(controlcard,'PDBREF').gt.0)
1333       iscode=index(controlcard,'ONE_LETTER')
1334       tree=(index(controlcard,'MAKE_TREE').gt.0)
1335       min_var=(index(controlcard,'MINVAR').gt.0)
1336       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1337       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1338       call readi(controlcard,'NCUT',ncut,1)
1339       call readi(controlcard,'SYM',symetr,1)
1340       write (iout,*) 'sym', symetr
1341       call readi(controlcard,'NSTART',nstart,0)
1342       call reada(controlcard,'BOXX',boxxsize,100.0d0)
1343       call reada(controlcard,'BOXY',boxysize,100.0d0)
1344       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1345       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
1346       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
1347       if (lipthick.gt.0.0d0) then
1348        bordliptop=(boxzsize+lipthick)/2.0
1349        bordlipbot=bordliptop-lipthick
1350       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
1351       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
1352       buflipbot=bordlipbot+lipbufthick
1353       bufliptop=bordliptop-lipbufthick
1354       if ((lipbufthick*2.0d0).gt.lipthick) &
1355        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
1356       endif !lipthick.gt.0
1357       write(iout,*) "bordliptop=",bordliptop
1358       write(iout,*) "bordlipbot=",bordlipbot
1359       write(iout,*) "bufliptop=",bufliptop
1360       write(iout,*) "buflipbot=",buflipbot
1361
1362       call readi(controlcard,'NCLUST',nclust,5)
1363 !      ions=index(controlcard,"IONS").gt.0
1364       call readi(controlcard,"SCELEMODE",scelemode,0)
1365       print *,"SCELE",scelemode
1366       call readi(controlcard,'TORMODE',tor_mode,0)
1367 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1368         write(iout,*) "torsional and valence angle mode",tor_mode
1369       call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
1370       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1371       write(iout,*) "R_CUT_ELE=",r_cut_ele
1372       call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
1373       call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
1374       call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
1375       call readi(controlcard,'NEND',nend,0)
1376       call reada(controlcard,'ECUT',ecut,10.0d0)
1377       call reada(controlcard,'PROB',prob_limit,0.99d0)
1378       write (iout,*) "Probability limit",prob_limit
1379       lgrp=(index(controlcard,'LGRP').gt.0)
1380       caonly=(index(controlcard,'CA_ONLY').gt.0)
1381       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1382       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1383       call readi(controlcard,'IOPT',iopt,2)
1384       lside = index(controlcard,"SIDE").gt.0
1385       efree = index(controlcard,"EFREE").gt.0
1386       call readi(controlcard,'NTEMP',nT,1)
1387       write (iout,*) "nT",nT
1388 !el      call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1389       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1390       write (iout,*) "nT",nT
1391       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1392       do i=1,nT
1393         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1394       enddo
1395       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1396       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1397       lprint_int=index(controlcard,"PRINT_INT") .gt.0
1398       if (min_var) iopt=1
1399       return
1400       end subroutine read_control
1401 !-----------------------------------------------------------------------------
1402       subroutine molread
1403 !
1404 ! Read molecular data.
1405 !
1406       use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
1407       use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1408 !                 wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1409 !                 wturn3,wturn4,wturn6,wvdwpp,weights
1410       use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1411                  indpdb
1412       use geometry, only: chainbuild,alloc_geo_arrays
1413       use energy, only: alloc_ener_arrays
1414       use io_base, only: rescode
1415       use control, only: setup_var,init_int_table
1416       use conform_compar, only: contact
1417 !      implicit none
1418 !      include 'DIMENSIONS'
1419 !      include 'COMMON.IOUNITS'
1420 !      include 'COMMON.GEO'
1421 !      include 'COMMON.VAR'
1422 !      include 'COMMON.INTERACT'
1423 !      include 'COMMON.LOCAL'
1424 !      include 'COMMON.NAMES'
1425 !      include 'COMMON.CHAIN'
1426 !      include 'COMMON.FFIELD'
1427 !      include 'COMMON.SBRIDGE'
1428 !      include 'COMMON.HEADER'
1429 !      include 'COMMON.CONTROL'
1430 !      include 'COMMON.CONTACTS'
1431 !      include 'COMMON.TIME1'
1432 !#ifdef MPL
1433 !      include 'COMMON.INFO'
1434 !#endif
1435       character(len=4) ::  sequence(nres) !(maxres)
1436       character(len=800) :: weightcard
1437 !      integer rescode
1438       real(kind=8) :: x(6*nres) !(maxvar)
1439       real(kind=8) :: ssscale
1440       integer  :: itype_pdb(nres) !(maxres)
1441 !      logical seq_comp
1442       integer :: i,j,kkk,mnum
1443 !
1444 ! Body
1445 !
1446 !el      allocate(weights(n_ene))
1447       allocate(weights(max_ene))
1448       call alloc_geo_arrays
1449       call alloc_ener_arrays
1450 !-----------------------------
1451       allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1452       allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1453       allocate(itype(nres+2,5)) !(maxres)
1454       allocate(molnum(nres+1))
1455       allocate(itel(nres+2))
1456
1457       do i=1,2*nres+2
1458         do j=1,3
1459           c(j,i)=0
1460           dc(j,i)=0
1461         enddo
1462       enddo
1463       do j=1,5
1464       do i=1,nres+2
1465         itype(i,j)=0
1466         itel(i)=0
1467       enddo
1468       enddo
1469 !--------------------------
1470 ! Read weights of the subsequent energy terms.
1471       call card_concat(weightcard,.true.)
1472       call reada(weightcard,'WSC',wsc,1.0d0)
1473       call reada(weightcard,'WLONG',wsc,wsc)
1474       call reada(weightcard,'WSCP',wscp,1.0d0)
1475       call reada(weightcard,'WELEC',welec,1.0D0)
1476       call reada(weightcard,'WVDWPP',wvdwpp,welec)
1477       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1478       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1479       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1480       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1481       call reada(weightcard,'WTURN3',wturn3,1.0D0)
1482       call reada(weightcard,'WTURN4',wturn4,1.0D0)
1483       call reada(weightcard,'WTURN6',wturn6,1.0D0)
1484       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1485       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1486       call reada(weightcard,'WBOND',wbond,1.0D0)
1487       call reada(weightcard,'WTOR',wtor,1.0D0)
1488       call reada(weightcard,'WTORD',wtor_d,1.0D0)
1489       call reada(weightcard,'WANG',wang,1.0D0)
1490       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1491       call reada(weightcard,'SCAL14',scal14,0.4D0)
1492       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1493       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1494       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1495       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
1496       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
1497       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
1498       if (index(weightcard,'SOFT').gt.0) ipot=6
1499       call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
1500       call reada(weightcard,'WELPP',welpp,0.0d0)
1501       call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
1502       call reada(weightcard,'WELPSB',welpsb,0.0D0)
1503       call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
1504       call reada(weightcard,'WELSB',welsb,0.0D0)
1505       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
1506       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
1507       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
1508       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
1509 !      print *,"KUR..",wtor_nucl
1510       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
1511       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
1512       call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
1513 ! 12/1/95 Added weight for the multi-body term WCORR
1514       call reada(weightcard,'WCORRH',wcorr,1.0D0)
1515       call reada(weightcard,'WSCBASE',wscbase,0.0D0)
1516       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
1517       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
1518       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
1519       call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
1520       call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
1521       call reada(weightcard,"D0CM",d0cm,3.78d0)
1522       call reada(weightcard,"AKCM",akcm,15.1d0)
1523       call reada(weightcard,"AKTH",akth,11.0d0)
1524       call reada(weightcard,"AKCT",akct,12.0d0)
1525       call reada(weightcard,"V1SS",v1ss,-1.08d0)
1526       call reada(weightcard,"V2SS",v2ss,7.61d0)
1527       call reada(weightcard,"V3SS",v3ss,13.7d0)
1528       call reada(weightcard,"EBR",ebr,-5.50D0)
1529       call reada(weightcard,"ATRISS",atriss,0.301D0)
1530       call reada(weightcard,"BTRISS",btriss,0.021D0)
1531       call reada(weightcard,"CTRISS",ctriss,1.001D0)
1532       call reada(weightcard,"DTRISS",dtriss,1.001D0)
1533       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
1534       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
1535       call reada(weightcard,"LIPSCALE",lipscale,1.0D0)
1536       call reada(weightcard,"WTL",wliptran,1.0D0)
1537
1538       call reada(weightcard,"HT",Ht,0.0D0)
1539       if (dyn_ss) then
1540        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
1541         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
1542         akcm=akcm/wsc*ssscale
1543         akth=akth/wsc*ssscale
1544         akct=akct/wsc*ssscale
1545         v1ss=v1ss/wsc*ssscale
1546         v2ss=v2ss/wsc*ssscale
1547         v3ss=v3ss/wsc*ssscale
1548       else
1549         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
1550       endif
1551
1552       if (wcorr4.gt.0.0d0) wcorr=wcorr4
1553       weights(1)=wsc
1554       weights(2)=wscp
1555       weights(3)=welec
1556       weights(4)=wcorr
1557       weights(5)=wcorr5
1558       weights(6)=wcorr6
1559       weights(7)=wel_loc
1560       weights(8)=wturn3
1561       weights(9)=wturn4
1562       weights(10)=wturn6
1563       weights(11)=wang
1564       weights(12)=wscloc
1565       weights(13)=wtor
1566       weights(14)=wtor_d
1567       weights(15)=wstrain
1568       weights(16)=wvdwpp
1569       weights(17)=wbond
1570       weights(18)=scal14
1571 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
1572       weights(21)=wsccor
1573           weights(26)=wvdwpp_nucl
1574           weights(27)=welpp
1575           weights(28)=wvdwpsb
1576           weights(29)=welpsb
1577           weights(30)=wvdwsb
1578           weights(31)=welsb
1579           weights(32)=wbond_nucl
1580           weights(33)=wang_nucl
1581           weights(34)=wsbloc
1582           weights(35)=wtor_nucl
1583           weights(36)=wtor_d_nucl
1584           weights(37)=wcorr_nucl
1585           weights(38)=wcorr3_nucl
1586           weights(41)=wcatcat
1587           weights(42)=wcatprot
1588           weights(46)=wscbase
1589           weights(47)=wscpho
1590           weights(48)=wpeppho
1591           weights(49)=wpeppho
1592           weights(50)=wcatnucl
1593           weights(56)=wcat_tran
1594           weights(22)=wliptran
1595
1596       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1597         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1598         wturn4,wturn6,wsccor
1599    10 format (/'Energy-term weights (unscaled):'// &
1600        'WSCC=   ',f10.6,' (SC-SC)'/ &
1601        'WSCP=   ',f10.6,' (SC-p)'/ &
1602        'WELEC=  ',f10.6,' (p-p electr)'/ &
1603        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1604        'WBOND=  ',f10.6,' (stretching)'/ &
1605        'WANG=   ',f10.6,' (bending)'/ &
1606        'WSCLOC= ',f10.6,' (SC local)'/ &
1607        'WTOR=   ',f10.6,' (torsional)'/ &
1608        'WTORD=  ',f10.6,' (double torsional)'/ &
1609        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1610        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1611        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1612        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1613        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1614        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1615        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1616        'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1617        'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1618
1619       if (wcorr4.gt.0.0d0) then
1620         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1621          'between contact pairs of peptide groups'
1622         write (iout,'(2(a,f5.3/))') &
1623         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1624         'Range of quenching the correlation terms:',2*delt_corr
1625       else if (wcorr.gt.0.0d0) then
1626         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1627          'between contact pairs of peptide groups'
1628       endif
1629       write (iout,'(a,f8.3)') &
1630         'Scaling factor of 1,4 SC-p interactions:',scal14
1631       write (iout,'(a,f8.3)') &
1632         'General scaling factor of SC-p interactions:',scalscp
1633       r0_corr=cutoff_corr-delt_corr
1634       do i=1,20
1635         aad(i,1)=scalscp*aad(i,1)
1636         aad(i,2)=scalscp*aad(i,2)
1637         bad(i,1)=scalscp*bad(i,1)
1638         bad(i,2)=scalscp*bad(i,2)
1639       enddo
1640
1641       call flush(iout)
1642       print *,'indpdb=',indpdb,' pdbref=',pdbref
1643
1644 ! Read sequence if not taken from the pdb file.
1645       if (iscode.gt.0) then
1646         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1647       else
1648         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1649       endif
1650       print *,nres_molec(:),nres
1651 ! Convert sequence to numeric code
1652       do i=1,nres_molec(1)
1653         mnum=1
1654         molnum(i)=1
1655         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1656       enddo
1657       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
1658         mnum=2
1659         molnum(i)=2
1660         write (iout,*),i,sequence(i)
1661         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1662
1663       enddo
1664       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
1665         mnum=5
1666         molnum(i)=5
1667         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1668       enddo
1669       print *,nres
1670       print '(20i4)',(itype(i,molnum(i)),i=1,nres)
1671
1672       do i=1,nres-1
1673 #ifdef PROCOR
1674         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
1675             .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
1676 #else
1677         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
1678 #endif
1679           itel(i)=0
1680 #ifdef PROCOR
1681         else if (iabs(itype(i+1,1)).ne.20) then
1682 #else
1683         else if (iabs(itype(i,1)).ne.20) then
1684 #endif
1685           itel(i)=1
1686         else
1687           itel(i)=2
1688         endif
1689       enddo
1690       write (iout,*) "ITEL"
1691       do i=1,nres-1
1692         write (iout,*) i,itype(i,molnum(i)),itel(i)
1693       enddo
1694
1695       print *,'Call Read_Bridge.'
1696       call read_bridge
1697       nnt=1
1698       nct=nres
1699       print *,'NNT=',NNT,' NCT=',NCT
1700       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1701       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1702       if (nstart.lt.nnt) nstart=nnt
1703       if (nend.gt.nct .or. nend.eq.0) nend=nct
1704       write (iout,*) "nstart",nstart," nend",nend
1705       nres0=nres
1706 !      if (pdbref) then
1707 !        read(inp,'(a)') pdbfile
1708 !        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1709 !        open(ipdbin,file=pdbfile,status='old',err=33)
1710 !        goto 34
1711 !  33    write (iout,'(a)') 'Error opening PDB file.'
1712 !        stop
1713 !  34    continue
1714 !        print *,'Begin reading pdb data'
1715 !        call readpdb
1716 !        print *,'Finished reading pdb data'
1717 !        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1718 !        do i=1,nres
1719 !          itype_pdb(i)=itype(i)
1720 !        enddo
1721 !        close (ipdbin)
1722 !        write (iout,'(a,i3)') 'nsup=',nsup
1723 !        nstart_seq=nnt
1724 !        if (nsup.le.(nct-nnt+1)) then
1725 !          do i=0,nct-nnt+1-nsup
1726 !            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1727 !              nstart_seq=nnt+i
1728 !              goto 111
1729 !            endif
1730 !          enddo
1731 !          write (iout,'(a)')
1732 !     &            'Error - sequences to be superposed do not match.'
1733 !          stop
1734 !        else
1735 !          do i=0,nsup-(nct-nnt+1)
1736 !            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1737 !     &      then
1738 !              nstart_sup=nstart_sup+i
1739 !              nsup=nct-nnt+1
1740 !              goto 111
1741 !            endif
1742 !          enddo
1743 !          write (iout,'(a)')
1744 !     &            'Error - sequences to be superposed do not match.'
1745 !        endif
1746 !  111   continue
1747 !        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1748 !     &                 ' nstart_seq=',nstart_seq
1749 !      endif
1750 write(iout,*)"przed ini_int_tab"
1751       call init_int_table
1752 write(iout,*)"po ini_int_tab"
1753 write(iout,*)"przed setup var"
1754       call setup_var
1755 write(iout,*)"po setup var"
1756       if (ns.gt.0.and.dyn_ss) then
1757           do i=nss+1,nhpb
1758             ihpb(i-nss)=ihpb(i)
1759             jhpb(i-nss)=jhpb(i)
1760             forcon(i-nss)=forcon(i)
1761             dhpb(i-nss)=dhpb(i)
1762           enddo
1763           nhpb=nhpb-nss
1764           nss=0
1765           link_start=1
1766           link_end=nhpb
1767 !          call hpb_partition
1768           do i=1,ns
1769             dyn_ss_mask(iss(i))=.true.
1770           enddo
1771       endif
1772
1773       write (iout,*) "molread: REFSTR",refstr
1774       if (refstr) then
1775         if (.not.pdbref) then
1776           call read_angles(inp,*38)
1777           goto 39
1778    38     write (iout,'(a)') 'Error reading reference structure.'
1779 #ifdef MPL
1780           call mp_stopall(Error_Msg)
1781 #else
1782           stop 'Error reading reference structure'
1783 #endif
1784    39     call chainbuild
1785           nstart_sup=nnt
1786           nstart_seq=nnt
1787           nsup=nct-nnt+1
1788           kkk=1
1789           do i=1,2*nres
1790             do j=1,3
1791               cref(j,i,kkk)=c(j,i)
1792             enddo
1793           enddo
1794         endif
1795         call contact(.true.,ncont_ref,icont_ref)
1796       endif
1797       return
1798       end subroutine molread
1799 !-----------------------------------------------------------------------------
1800       subroutine openunits
1801 !      implicit none
1802 !      include 'DIMENSIONS'
1803       use control_data, only: from_cx,from_bx,from_cart
1804 #ifdef MPI
1805       use MPI_data
1806       include "mpif.h"
1807       character(len=3) :: liczba
1808 !      include "COMMON.MPI"
1809 #endif
1810 !      include 'COMMON.IOUNITS'
1811 !      include 'COMMON.CONTROL'
1812       integer :: lenpre,lenpot !,ilen
1813 !      external ilen
1814       character(len=16) :: cformat,cprint
1815 !      character(len=16) ucase
1816       integer :: lenint,lenout
1817       call getenv('INPUT',prefix)
1818       call getenv('OUTPUT',prefout)
1819       call getenv('INTIN',prefintin)
1820       call getenv('COORD',cformat)
1821       call getenv('PRINTCOOR',cprint)
1822       call getenv('SCRATCHDIR',scratchdir)
1823       from_bx=.true.
1824       from_cx=.false.
1825       if (index(ucase(cformat),'CX').gt.0) then
1826         from_cx=.true.
1827         from_bx=.false.
1828       endif
1829       from_cart=.true.
1830       lenpre=ilen(prefix)
1831       lenout=ilen(prefout)
1832       lenint=ilen(prefintin)
1833 ! Get the names and open the input files
1834       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1835 #ifdef MPI
1836       write (liczba,'(bz,i3.3)') me
1837       outname=prefout(:lenout)//'_clust.out_'//liczba
1838 #else
1839       outname=prefout(:lenout)//'_clust.out'
1840 #endif
1841       if (from_bx) then
1842         intinname=prefintin(:lenint)//'.bx'
1843       else if (from_cx) then
1844         intinname=prefintin(:lenint)//'.cx'
1845       else
1846         intinname=prefintin(:lenint)//'.int'
1847       endif
1848       rmsname=prefintin(:lenint)//'.rms'
1849       open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1850              status='unknown')
1851       open (jrms,file=rmsname,status='unknown')
1852       open(iout,file=outname,status='unknown')
1853 ! Get parameter filenames and open the parameter files.
1854       call getenv('BONDPAR',bondname)
1855       open (ibond,file=bondname,status='old')
1856       call getenv('THETPAR',thetname)
1857       open (ithep,file=thetname,status='old')
1858       call getenv('ROTPAR',rotname)
1859       open (irotam,file=rotname,status='old')
1860       call getenv('TORPAR',torname)
1861       open (itorp,file=torname,status='old')
1862       call getenv('TORDPAR',tordname)
1863       open (itordp,file=tordname,status='old')
1864       call getenv('FOURIER',fouriername)
1865       open (ifourier,file=fouriername,status='old')
1866       call getenv('ELEPAR',elename)
1867       open (ielep,file=elename,status='old')
1868       call getenv('SIDEPAR',sidename)
1869       open (isidep,file=sidename,status='old')
1870       call getenv('SIDEP',sidepname)
1871       open (isidep1,file=sidepname,status="old")
1872       call getenv('SCCORPAR',sccorname)
1873       open (isccor,file=sccorname,status="old")
1874       call getenv('BONDPAR_NUCL',bondname_nucl)
1875       open (ibond_nucl,file=bondname_nucl,status='old')
1876       call getenv('THETPAR_NUCL',thetname_nucl)
1877       open (ithep_nucl,file=thetname_nucl,status='old')
1878       call getenv('ROTPAR_NUCL',rotname_nucl)
1879       open (irotam_nucl,file=rotname_nucl,status='old')
1880       call getenv('TORPAR_NUCL',torname_nucl)
1881       open (itorp_nucl,file=torname_nucl,status='old')
1882       call getenv('TORDPAR_NUCL',tordname_nucl)
1883       open (itordp_nucl,file=tordname_nucl,status='old')
1884       call getenv('SIDEPAR_NUCL',sidename_nucl)
1885       open (isidep_nucl,file=sidename_nucl,status='old')
1886       call getenv('SIDEPAR_SCBASE',sidename_scbase)
1887       open (isidep_scbase,file=sidename_scbase,status='old')
1888       call getenv('PEPPAR_PEPBASE',pepname_pepbase)
1889       open (isidep_pepbase,file=pepname_pepbase,status='old')
1890       call getenv('SCPAR_PHOSPH',pepname_scpho)
1891       open (isidep_scpho,file=pepname_scpho,status='old')
1892       call getenv('PEPPAR_PHOSPH',pepname_peppho)
1893       open (isidep_peppho,file=pepname_peppho,status='old')
1894
1895
1896       call getenv('LIPTRANPAR',liptranname)
1897       open (iliptranpar,file=liptranname,status='old')
1898       call getenv('TUBEPAR',tubename)
1899       open (itube,file=tubename,status='old')
1900       call getenv('IONPAR',ionname)
1901       open (iion,file=ionname,status='old')
1902       call getenv('IONPAR_NUCL',ionnuclname)
1903       open (iionnucl,file=ionnuclname,status='old')
1904       call getenv('IONPAR_TRAN',iontranname)
1905       open (iiontran,file=iontranname,status='old')
1906
1907
1908 #ifndef OLDSCP
1909 !
1910 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1911 ! Use -DOLDSCP to use hard-coded constants instead.
1912 !
1913       call getenv('SCPPAR',scpname)
1914       open (iscpp,file=scpname,status='old')
1915       call getenv('SCPPAR_NUCL',scpname_nucl)
1916       open (iscpp_nucl,file=scpname_nucl,status='old')
1917
1918 #endif
1919       return
1920       end subroutine openunits
1921 !-----------------------------------------------------------------------------
1922 ! geomout.F
1923 !-----------------------------------------------------------------------------
1924       subroutine pdboutC(etot,rmsd,tytul)
1925
1926       use energy_data, only: ihpb,jhpb,itype,molnum
1927       use energy, only:boxshift,to_box
1928 !      implicit real*8 (a-h,o-z)
1929 !      include 'DIMENSIONS'
1930 !      include 'COMMON.CONTROL'
1931 !      include 'COMMON.CHAIN'
1932 !      include 'COMMON.INTERACT'
1933 !      include 'COMMON.NAMES'
1934 !      include 'COMMON.IOUNITS'
1935 !      include 'COMMON.HEADER'
1936 !      include 'COMMON.SBRIDGE'
1937 !      include 'COMMON.TEMPFAC'
1938       character(len=50) :: tytul
1939       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1940                                                'G','H','I','J'/)
1941       integer :: ica(nres),k,iti1,ki
1942       real(kind=8) :: etot,rmsd
1943       integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
1944       real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
1945       write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1946         ' ENERGY ',etot,' RMS ',rmsd
1947       iatom=0
1948       ichain=1
1949       ires=0
1950 !      boxxxshift(1)=int(c(1,nnt)/boxxsize)
1951 !      if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
1952 !      boxxxshift(2)=int(c(2,nnt)/boxzsize)
1953 !      if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
1954 !      boxxxshift(3)=int(c(3,nnt)/boxzsize)
1955 !      if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
1956 !      do i=1,3
1957 !        if (boxxxshift(i).gt.2) boxxxshift(i)=2
1958 !        if (boxxxshift(i).lt.-2) boxxxshift(i)=-2
1959 !
1960 !      enddo
1961       do j=1,3
1962       cbeg(j)=c(j,nnt)
1963       enddo
1964       call to_box(cbeg(1),cbeg(2),cbeg(3))
1965       boxxxx(1)=boxxsize
1966       boxxxx(2)=boxysize
1967       boxxxx(3)=boxzsize
1968
1969       do i=nnt,nct
1970         mnum=molnum(i)
1971         iti=itype(i,mnum)
1972         if (i.ne.nct) then
1973         iti1=itype(i+1,mnum)
1974         mnum1=molnum(i+1)
1975         if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle
1976         endif
1977         if (iti.eq.ntyp1_molec(mnum)) then
1978           ichain=ichain+1
1979           ires=0
1980           write (ipdb,'(a)') 'TER'
1981         else
1982         ires=ires+1
1983 !         if (molnum(i).ge.1) then
1984           if (i.le.3) then
1985            ki=2
1986           else
1987           if (itype(i,mnum).ne.ntyp1_molec(mnum)) then
1988           ki=ki+1
1989           endif
1990          endif
1991 !        endif
1992 !          print *,"tu2",i,k
1993           if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
1994           if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
1995         do j=1,3
1996         cbeg(j)=c(j,ki)
1997         enddo
1998         iatom=iatom+1
1999         ica(i)=iatom
2000         call to_box(c(1,i),c(2,i),c(3,i))
2001
2002         DO k=1,3
2003         Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k))
2004         c(k,i)=cbeg(k)+Rdist(k)
2005         ENDDO
2006       if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then
2007         call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres))
2008         DO k=1,3
2009         Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k))
2010         c(k,i+nres)=cbeg(k)+Rdist(k)
2011         ENDDO
2012        endif
2013         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
2014            ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
2015         if ((iti.ne.10).and.(mnum.ne.5)) then
2016           iatom=iatom+1
2017           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
2018             ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
2019         endif
2020         endif
2021       enddo
2022       write (ipdb,'(a)') 'TER'
2023       do i=nnt,nct-1
2024         mnum=molnum(i)
2025         mnum1=molnum(i+1)
2026         if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
2027         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2028           write (ipdb,30) ica(i),ica(i+1)
2029         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2030           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2031         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
2032           write (ipdb,30) ica(i),ica(i)+1
2033         endif
2034       enddo
2035       if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
2036         write (ipdb,30) ica(nct),ica(nct)+1
2037       endif
2038       do i=1,nss
2039         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2040       enddo
2041       write (ipdb,'(a6)') 'ENDMDL'
2042   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2043   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2044   30  FORMAT ('CONECT',8I5)
2045       return
2046       end subroutine pdboutC
2047 !-----------------------------------------------------------------------------
2048       subroutine cartout(igr,i,etot,free,rmsd,plik)
2049 !     implicit real*8 (a-h,o-z)
2050 !     include 'DIMENSIONS'
2051 !     include 'sizesclu.dat'
2052 !     include 'COMMON.IOUNITS'
2053 !     include 'COMMON.CHAIN'
2054 !     include 'COMMON.VAR'
2055 !     include 'COMMON.LOCAL'
2056 !     include 'COMMON.INTERACT'
2057 !     include 'COMMON.NAMES'
2058 !     include 'COMMON.GEO'
2059 !     include 'COMMON.CLUSTER'
2060       integer :: igr,i,j,k
2061       real(kind=8) :: etot,free,rmsd
2062       character(len=80) :: plik
2063       open (igeom,file=plik,position='append')
2064       write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
2065       write (igeom,'(i4,$)') &
2066         nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
2067       write (igeom,'(i10)') iscore(i)
2068       write (igeom,'(8f10.5)') &
2069         ((allcart(k,j,i),k=1,3),j=1,nres),&
2070         ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
2071       return
2072       end subroutine cartout
2073 !------------------------------------------------------------------------------
2074 !      subroutine alloc_clust_arrays(n_conf)
2075
2076 !      integer :: n_conf
2077 !COMMON.CLUSTER
2078 !      common /clu/
2079 !      allocate(diss(maxdist)) !(maxdist)
2080 !el      allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
2081 !      allocatable :: enetb !(max_ene,maxstr_proc)
2082 !el      allocate(entfac(maxconf)) !(maxconf)
2083 !      allocatable :: totfree_gr !(maxgr)
2084 !el      allocate(rcutoff(max_cut+1)) !(max_cut+1)
2085 !      common /clu1/
2086 !      allocatable :: licz,iass !(maxgr)
2087 !      allocatable :: nconf !(maxgr,maxingr)
2088 !      allocatable :: iass_tot !(maxgr,max_cut)
2089 !      allocatable :: list_conf !(maxconf)
2090 !      common /alles/
2091 !el      allocatable :: allcart !(3,maxres2,maxstr_proc)
2092 !el      allocate(rmstb(maxconf)) !(maxconf)
2093 !el      allocate(mult(nres)) !(maxres)
2094 !el      allocatable :: nss_all !(maxstr_proc)
2095 !el      allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
2096 !      allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
2097 !COMMON.TEMPFAC
2098 !      common /factemp/
2099 !      allocatable :: tempfac !(2,maxres)
2100 !COMMON.FREE
2101 !      common /free/
2102 !el      allocate(beta_h(maxT)) !(maxT)
2103
2104 !      end subroutine alloc_clust_arrays
2105 !-----------------------------------------------------------------------------
2106 !-----------------------------------------------------------------------------
2107       end module io_clust