source/cluster/io_clust.F90

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