--- /dev/null
+ module io_clust
+!-----------------------------------------------------------------------------
+ use clust_data
+ use io_units
+! use names
+ use io_base !, only: ilen
+ use geometry_data, only: nres,c
+ use energy_data, only: nnt,nct,nss
+ use control_data, only: lside
+ implicit none
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+ contains
+!-----------------------------------------------------------------------------
+! wrtclust.f
+!-----------------------------------------------------------------------------
+ SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
+
+ use hc_, only: ioffset
+ use control_data, only: lprint_cart,lprint_int,titel
+ use geometry, only: int_from_cart1,nres
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'sizesclu.dat'
+ integer,parameter :: num_in_line=5
+ LOGICAL :: PRINTANG(max_cut)
+ integer :: PRINTPDB(max_cut),printmol2(max_cut)
+! include 'COMMON.CONTROL'
+! include 'COMMON.HEADER'
+! include 'COMMON.CHAIN'
+! include 'COMMON.VAR'
+! include 'COMMON.CLUSTER'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.FREE'
+! include 'COMMON.TEMPFAC'
+ real(kind=8) :: rmsave(maxgr)
+ CHARACTER(len=64) :: prefixp,NUMM,MUMM,EXTEN,extmol
+ character(len=80) :: cfname
+ character(len=8) :: ctemper
+ DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,&
+ MUMM /'000'/
+! external ilen
+ integer :: ncon,icut,ib
+ integer :: i,ii,ii1,ii2,igr,ind1,ind2,ico,icon,&
+ irecord,nrecord,j,k,jj,ind,ncon_lim,ncon_out
+ real(kind=8) :: temper,curr_dist,emin,qpart,boltz,&
+ ave_dim,amax_dim,emin1
+
+
+ allocate(tempfac(2,nres))
+
+ do i=1,64
+ cfname(i:i)=" "
+ enddo
+! print *,"calling WRTCLUST",ncon
+! write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
+ rewind 80
+ call flush(iout)
+ temper=1.0d0/(beta_h(ib)*1.987d-3)
+ if (temper.lt.100.0d0) then
+ write(ctemper,'(f3.0)') temper
+ ctemper(3:3)=" "
+ else if (temper.lt.1000.0) then
+ write (ctemper,'(f4.0)') temper
+ ctemper(4:4)=" "
+ else
+ write (ctemper,'(f5.0)') temper
+ ctemper(5:5)=" "
+ endif
+
+ do i=1,ncon*(ncon-1)/2
+ read (80) diss(i)
+ enddo
+ close(80,status='delete')
+!
+! PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
+!
+ ii1= index(intinname,'/')
+ ii2=ii1
+ ii1=ii1+1
+ do while (ii2.gt.0)
+ ii1=ii1+ii2
+ ii2=index(intinname(ii1:),'/')
+ enddo
+ ii = ii1+index(intinname(ii1:),'.')-1
+ if (ii.eq.0) then
+ ii=ilen(intinname)
+ else
+ ii=ii-1
+ endif
+ prefixp=intinname(ii1:ii)
+!d print *,icut,printang(icut),printpdb(icut),printmol2(icut)
+!d print *,'ecut=',ecut
+ WRITE (iout,100) NGR
+ DO 19 IGR=1,NGR
+ WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
+ NRECORD=LICZ(IGR)/num_in_line
+ IND1=1
+ DO 63 IRECORD=1,NRECORD
+ IND2=IND1+num_in_line-1
+ WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
+ totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
+ IND1=IND2+1
+ 63 CONTINUE
+ WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
+ totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
+ IND1=1
+ ICON=list_conf(NCONF(IGR,1))
+! WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
+! 12/8/93 Estimation of "diameters" of the subsequent families.
+ ave_dim=0.0
+ amax_dim=0.0
+! write (iout,*) "ecut",ecut
+ do i=2,licz(igr)
+ ii=nconf(igr,i)
+ if (totfree(ii)-emin .gt. ecut) goto 10
+ do j=1,i-1
+ jj=nconf(igr,j)
+ if (jj.eq.1) exit
+ if (ii.lt.jj) then
+ ind=ioffset(ncon,ii,jj)
+ else
+ ind=ioffset(ncon,jj,ii)
+ endif
+! write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
+! & " ind",ind
+ call flush(iout)
+ curr_dist=dabs(diss(ind)+0.0d0)
+! write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
+! & list_conf(jj),curr_dist
+ if (curr_dist .gt. amax_dim) amax_dim=curr_dist
+ ave_dim=ave_dim+curr_dist**2
+ enddo
+ enddo
+ 10 if (licz(igr) .gt. 1) &
+ ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
+ write (iout,'(/A,F8.1,A,F8.1)') &
+ 'Max. distance in the family:',amax_dim,&
+ '; average distance in the family:',ave_dim
+ rmsave(igr)=0.0d0
+ qpart=0.0d0
+ do i=1,licz(igr)
+ icon=nconf(igr,i)
+ boltz=dexp(-totfree(icon))
+ rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
+ qpart=qpart+boltz
+ enddo
+ rmsave(igr)=rmsave(igr)/qpart
+ write (iout,'(a,f5.2,a)') "Average RMSD",rmsave(igr)," A"
+ 19 CONTINUE
+ WRITE (iout,400)
+ WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
+! print *,icut,printang(icut)
+ IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
+ emin=totfree_gr(1)
+! print *,'emin',emin,' ngr',ngr
+ if (lprint_cart) then
+ cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
+ //"K"//".x"
+ else
+ cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
+ //"K"//".int"
+ endif
+ do igr=1,ngr
+ icon=nconf(igr,1)
+ if (totfree_gr(igr)-emin.le.ecut) then
+ if (lprint_cart) then
+ call cartout(igr,icon,totfree(icon)/beta_h(ib),&
+ totfree_gr(igr)/beta_h(ib),&
+ rmstb(icon),cfname)
+ else
+! print '(a)','calling briefout'
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=allcart(j,i,icon)
+ enddo
+ enddo
+ call int_from_cart1(.false.)
+!el call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),&
+!el totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),&
+!el jhpb_all(1,icon),cfname)
+ call briefout(igr,totfree(icon)/beta_h(ib),&
+ totfree_gr(igr))
+! print '(a)','exit briefout'
+ endif
+ endif
+ enddo
+ close(igeom)
+ ENDIF
+ IF (PRINTPDB(ICUT).gt.0) THEN
+! Write out a number of conformations from each family in PDB format and
+! create InsightII command file for their displaying in different colors
+ cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
+ //"K_"//'ave'//exten
+ write (iout,*) "cfname",cfname
+ OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
+ write (ipdb,'(a,f8.2)') &
+ "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
+ close (ipdb)
+ I=1
+ ICON=NCONF(1,1)
+ EMIN=totfree_gr(I)
+ emin1=totfree(icon)
+ DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
+! write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
+ write (NUMM,'(bz,i4.4)') i
+ ncon_lim=min0(licz(i),printpdb(icut))
+ cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
+ //"K_"//numm(:ilen(numm))//exten
+ OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
+ write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5," AVE RMSD",0pf5.2)')&
+ i,totfree_gr(i)/beta_h(ib),rmsave(i) !'
+! Write conformations of the family i to PDB files
+ ncon_out=1
+ do while (ncon_out.lt.printpdb(icut) .and. &
+ ncon_out.lt.licz(i).and. &
+ totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
+ ncon_out=ncon_out+1
+! write (iout,*) i,ncon_out,nconf(i,ncon_out),
+! & totfree(nconf(i,ncon_out)),emin1,ecut
+ enddo
+ write (iout,*) "ncon_out",ncon_out
+ call flush(iout)
+ do j=1,nres
+ tempfac(1,j)=5.0d0
+ tempfac(2,j)=5.0d0
+ enddo
+ do j=1,ncon_out
+ icon=nconf(i,j)
+ do ii=1,2*nres
+ do k=1,3
+ c(k,ii)=allcart(k,ii,icon)
+ enddo
+ enddo
+ call pdboutC(totfree(icon)/beta_h(ib),rmstb(icon),titel)
+ write (ipdb,'("TER")')
+ enddo
+ close(ipdb)
+! Average structures and structures closest to average
+ cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
+ //"K_"//'ave'//exten
+ OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',&
+ position="APPEND")
+ call ave_coord(i)
+ write (ipdb,'(a,i5)') "REMARK CLUSTER",i
+ call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
+ write (ipdb,'("TER")')
+ call closest_coord(i)
+ call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
+ write (ipdb,'("TER")')
+ close (ipdb)
+ I=I+1
+ ICON=NCONF(I,1)
+ emin1=totfree(icon)
+ ENDDO
+ ENDIF
+ IF (printmol2(icut).gt.0) THEN
+! Write out a number of conformations from each family in PDB format and
+! create InsightII command file for their displaying in different colors
+ I=1
+ ICON=NCONF(1,1)
+ EMIN=ENERGY(ICON)
+ emin1=emin
+ DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
+ write (NUMM,'(bz,i4.4)') i
+ cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
+ //"K_"//numm(:ilen(numm))//extmol
+ OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
+ ncon_out=1
+ do while (ncon_out.lt.printmol2(icut) .and. &
+ ncon_out.lt.licz(i).and. &
+ totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
+ ncon_out=ncon_out+1
+ enddo
+ do j=1,ncon_out
+ icon=nconf(i,j)
+ do ii=1,2*nres
+ do k=1,3
+ c(k,ii)=allcart(k,ii,icon)
+ enddo
+ enddo
+ CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
+ enddo
+ CLOSE(imol2)
+ I=I+1
+ ICON=NCONF(I,1)
+ emin1=totfree(icon)
+ ENDDO
+ ENDIF
+ 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
+ 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,&
+ ' CONTAINS ',I4,' CONFORMATION(S): ')
+! 300 FORMAT ( 8(I4,F6.1))
+ 300 FORMAT (5(I4,1pe12.3))
+ 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
+ 500 FORMAT (8(2I4,2X))
+ 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
+ RETURN
+ END SUBROUTINE WRTCLUST
+!------------------------------------------------------------------------------
+ subroutine ave_coord(igr)
+
+ use control_data, only:lside
+ use regularize_, only:fitsq,matvec
+! implicit none
+! include 'DIMENSIONS'
+! include 'sizesclu.dat'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CLUSTER'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.TEMPFAC'
+! include 'COMMON.IOUNITS'
+ logical :: non_conv
+ real(kind=8) :: przes(3),obrot(3,3)
+ real(kind=8) :: xx(3,2*nres),yy(3,2*nres),csq(3,2*nres) !(3,maxres2)
+ real(kind=8) :: eref
+ integer :: i,ii,j,k,icon,jcon,igr
+ real(kind=8) :: rms,boltz,qpart,cwork(3,2*nres),cref1(3,2*nres) !(3,maxres2)
+! write (iout,*) "AVE_COORD: igr",igr
+ jcon=nconf(igr,1)
+ eref=totfree(jcon)
+ boltz = dexp(-totfree(jcon)+eref)
+ qpart=boltz
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=allcart(j,i,jcon)*boltz
+ cref1(j,i)=allcart(j,i,jcon)
+ csq(j,i)=allcart(j,i,jcon)**2*boltz
+ enddo
+ enddo
+ DO K=2,LICZ(IGR)
+ jcon=nconf(igr,k)
+ if (lside) then
+ ii=0
+ do i=nnt,nct
+ ii=ii+1
+ do j=1,3
+ xx(j,ii)=allcart(j,i,jcon)
+ yy(j,ii)=cref1(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+! if (itype(i).ne.10) then
+ ii=ii+1
+ do j=1,3
+ xx(j,ii)=allcart(j,i+nres,jcon)
+ yy(j,ii)=cref1(j,i+nres)
+ enddo
+! endif
+ enddo
+ call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
+ else
+ do i=nnt,nct
+ do j=1,3
+ cwork(j,i)=allcart(j,i,jcon)
+ enddo
+ enddo
+ call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot &
+ ,non_conv)
+ endif
+! write (iout,*) "rms",rms
+! do i=1,3
+! write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
+! enddo
+ if (rms.lt.0.0) then
+ print *,'error, rms^2 = ',rms,icon,jcon
+ stop
+ endif
+ if (non_conv) print *,non_conv,icon,jcon
+ boltz=dexp(-totfree(jcon)+eref)
+ qpart = qpart + boltz
+ do i=1,2*nres
+ do j=1,3
+ xx(j,i)=allcart(j,i,jcon)
+ enddo
+ enddo
+ call matvec(cwork,obrot,xx,2*nres)
+ do i=1,2*nres
+! write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
+! & (allcart(j,i,jcon),j=1,3)
+ do j=1,3
+ cwork(j,i)=cwork(j,i)+przes(j)
+ c(j,i)=c(j,i)+cwork(j,i)*boltz
+ csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
+ enddo
+ enddo
+ ENDDO ! K
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=c(j,i)/qpart
+ csq(j,i)=csq(j,i)/qpart-c(j,i)**2
+ enddo
+! write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
+ enddo
+ do i=nnt,nct
+ tempfac(1,i)=0.0d0
+ tempfac(2,i)=0.0d0
+ do j=1,3
+ tempfac(1,i)=tempfac(1,i)+csq(j,i)
+ tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
+ enddo
+ tempfac(1,i)=dsqrt(tempfac(1,i))
+ tempfac(2,i)=dsqrt(tempfac(2,i))
+ enddo
+ return
+ end subroutine ave_coord
+!------------------------------------------------------------------------------
+ subroutine closest_coord(igr)
+
+ use regularize_, only:fitsq
+! implicit none
+! include 'DIMENSIONS'
+! include 'sizesclu.dat'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CLUSTER'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+ logical :: non_conv
+ real(kind=8) :: przes(3),obrot(3,3)
+ real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
+ integer :: i,ii,j,k,icon,jcon,jconmin,igr
+ real(kind=8) :: rms,rmsmin,cwork(3,2*nres)
+ rmsmin=1.0d10
+ jconmin=nconf(igr,1)
+ DO K=1,LICZ(IGR)
+ jcon=nconf(igr,k)
+ if (lside) then
+ ii=0
+ do i=nnt,nct
+ ii=ii+1
+ do j=1,3
+ xx(j,ii)=allcart(j,i,jcon)
+ yy(j,ii)=c(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+! if (itype(i).ne.10) then
+ ii=ii+1
+ do j=1,3
+ xx(j,ii)=allcart(j,i+nres,jcon)
+ yy(j,ii)=c(j,i+nres)
+ enddo
+! endif
+ enddo
+ call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
+ else
+ do i=nnt,nct
+ do j=1,3
+ cwork(j,i)=allcart(j,i,jcon)
+ enddo
+ enddo
+ call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot,&
+ non_conv)
+ endif
+ if (rms.lt.0.0) then
+ print *,'error, rms^2 = ',rms,icon,jcon
+ stop
+ endif
+! write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
+ if (non_conv) print *,non_conv,icon,jcon
+ if (rms.lt.rmsmin) then
+ rmsmin=rms
+ jconmin=jcon
+ endif
+ ENDDO ! K
+! write (iout,*) "rmsmin",rmsmin," rms",rms
+ call flush(iout)
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=allcart(j,i,jconmin)
+ enddo
+ enddo
+ return
+ end subroutine closest_coord
+!-----------------------------------------------------------------------------
+! read_coords.F
+!-----------------------------------------------------------------------------
+ subroutine read_coords(ncon,*)
+
+ use energy_data, only: ihpb,jhpb,max_ene
+ use control_data, only: from_bx,from_cx
+ use control, only: tcpu
+! implicit none
+! include "DIMENSIONS"
+! include "sizesclu.dat"
+#ifdef MPI
+ use MPI_data
+ include "mpif.h"
+ integer :: IERROR,ERRCODE !,STATUS(MPI_STATUS_SIZE)
+! include "COMMON.MPI"
+#else
+ use MPI_data, only: nprocs
+#endif
+! include "COMMON.CONTROL"
+! include "COMMON.CHAIN"
+! include "COMMON.INTERACT"
+! include "COMMON.IOUNITS"
+! include "COMMON.VAR"
+! include "COMMON.SBRIDGE"
+! include "COMMON.GEO"
+! include "COMMON.CLUSTER"
+ character(len=3) :: liczba
+ integer :: ncon
+ integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,if,ib,&
+ nn,nn1,inan
+ integer :: ixdrf,iret,itmp
+ real(kind=4) :: prec,reini,refree,rmsdev
+ integer :: nrec,nlines,iscor,lenrec,lenrec_in
+ real(kind=8) :: energ,t_acq !,tcpu
+!el integer ilen,iroof
+!el external ilen,iroof
+ real(kind=8) :: rjunk
+ integer :: ntot_all(0:nprocs-1) !(0:maxprocs-1)
+ logical :: lerr
+ real(kind=8) :: energia(0:max_ene),etot
+ real(kind=4) :: csingle(3,2*nres+2)
+ integer :: Previous,Next
+ character(len=256) :: bprotfiles
+! print *,"Processor",me," calls read_protein_data"
+#ifdef MPI
+ if (me.eq.master) then
+ Previous=MPI_PROC_NULL
+ else
+ Previous=me-1
+ endif
+ if (me.eq.nprocs-1) then
+ Next=MPI_PROC_NULL
+ else
+ Next=me+1
+ endif
+! Set the scratchfile names
+ write (liczba,'(bz,i3.3)') me
+
+ allocate(STATUS(MPI_STATUS_SIZE))
+#endif
+! 1/27/05 AL Change stored coordinates to single precision and don't store
+! energy components in the binary databases.
+ lenrec=12*(nres+nct-nnt+1)+4*(2*nss+2)+16
+ lenrec_in=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
+#ifdef DEBUG
+ write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss", nss
+ write (iout,*) "lenrec_in",lenrec_in
+#endif
+ bprotfiles=scratchdir(:ilen(scratchdir))// &
+ "/"//prefix(:ilen(prefix))//liczba//".xbin"
+! EL
+! allocate cluster arrays
+ allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
+ allocate(entfac(maxconf)) !(maxconf)
+ allocate(rmstb(maxconf)) !(maxconf)
+ allocate(allcart(3,2*nres,maxstr_proc)) !(3,maxres2,maxstr_proc)
+ allocate(nss_all(maxstr_proc)) !(maxstr_proc)
+ allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc)
+ allocate(iscore(maxconf)) !(maxconf)
+
+
+#ifdef CHUJ
+ ICON=1
+ 123 continue
+ if (from_cart .and. .not. from_bx .and. .not. from_cx) then
+ if (efree) then
+ read (intin,*,end=13,err=11) energy(icon),totfree(icon),&
+ rmstb(icon),&
+ nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
+ i=1,nss_all(icon)),iscore(icon)
+ else
+ read (intin,*,end=13,err=11) energy(icon),rmstb(icon),&
+ nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
+ i=1,nss_all(icon)),iscore(icon)
+ endif
+ read (intin,'(8f10.5)',end=13,err=10) &
+ ((allcart(j,i,icon),j=1,3),i=1,nres),&
+ ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
+ print *,icon,energy(icon),nss_all(icon),rmstb(icon)
+ else
+ read(intin,'(a80)',end=13,err=12) lineh
+ read(lineh(:5),*,err=8) ic
+ if (efree) then
+ read(lineh(6:),*,err=8) energy(icon)
+ else
+ read(lineh(6:),*,err=8) energy(icon)
+ endif
+ goto 9
+ 8 ic=1
+ print *,'error, assuming e=1d10',lineh
+ energy(icon)=1d10
+ nss=0
+ 9 continue
+!old read(lineh(18:),*,end=13,err=11) nss_all(icon)
+ ii = index(lineh(15:)," ")+15
+ read(lineh(ii:),*,end=13,err=11) nss_all(icon)
+ IF (NSS_all(icon).LT.9) THEN
+ read (lineh(20:),*,end=102) &
+ (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),&
+ iscore(icon)
+ ELSE
+ read (lineh(20:),*,end=102) &
+ (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
+ read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),&
+ I=9,NSS_all(icon)),iscore(icon)
+ ENDIF
+
+ 102 continue
+
+ PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
+ call read_angles(intin,*13)
+ do i=1,nres
+ phiall(i,icon)=phi(i)
+ thetall(i,icon)=theta(i)
+ alphall(i,icon)=alph(i)
+ omall(i,icon)=omeg(i)
+ enddo
+ endif
+ ICON=ICON+1
+ GOTO 123
+!
+! CALCULATE DISTANCES
+!
+ 10 print *,'something wrong with angles'
+ goto 13
+ 11 print *,'something wrong with NSS',nss
+ goto 13
+ 12 print *,'something wrong with header'
+
+ 13 NCON=ICON-1
+
+#endif
+ call flush(iout)
+ jj_old=1
+ open (icbase,file=bprotfiles,status="unknown",&
+ form="unformatted",access="direct",recl=lenrec)
+! Read conformations from binary DA files (one per batch) and write them to
+! a binary DA scratchfile.
+ jj=0
+ jjj=0
+#ifdef MPI
+ write (liczba,'(bz,i3.3)') me
+ IF (ME.EQ.MASTER) THEN
+! Only the master reads the database; it'll send it to the other procs
+! through a ring.
+#endif
+ t_acq = tcpu()
+ icount=0
+
+ if (from_bx) then
+
+ open (intin,file=intinname,status="old",form="unformatted",&
+ access="direct",recl=lenrec_in)
+
+ else if (from_cx) then
+#if (defined(AIX) && !defined(JUBL))
+ call xdrfopen_(ixdrf,intinname, "r", iret)
+#else
+ call xdrfopen(ixdrf,intinname, "r", iret)
+#endif
+ prec=10000.0
+ write (iout,*) "xdrfopen: iret",iret
+ if (iret.eq.0) then
+ write (iout,*) "Error: coordinate file ",&
+ intinname(:ilen(intinname))," does not exist."
+ call flush(iout)
+#ifdef MPI
+ call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
+#endif
+ stop
+ endif
+ else
+ write (iout,*) "Error: coordinate format not specified"
+ call flush(iout)
+#ifdef MPI
+ call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
+#else
+ stop
+#endif
+ endif
+
+!#define DEBUG
+#ifdef DEBUG
+ write (iout,*) "Opening file ",intinname(:ilen(intinname))
+ write (iout,*) "lenrec",lenrec_in
+ call flush(iout)
+#endif
+!#undef DEBUG
+! write (iout,*) "maxconf",maxconf
+ i=0
+ do while (.true.)
+ i=i+1
+!el if (i.gt.maxconf) then
+!el write (iout,*) "Error: too many conformations ",&
+!el "(",maxconf,") maximum."
+!#ifdef MPI
+!el call MPI_Abort(MPI_COMM_WORLD,errcode,ierror)
+!#endif
+!el stop
+!el endif
+! write (iout,*) "i",i
+! call flush(iout)
+ if (from_bx) then
+ read(intin,err=101,end=101) &
+ ((csingle(l,k),l=1,3),k=1,nres),&
+ ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
+ nss,(ihpb(k),jhpb(k),k=1,nss),&
+ energy(jj+1),&
+ entfac(jj+1),rmstb(jj+1),iscor
+ do j=1,2*nres
+ do k=1,3
+ c(k,j)=csingle(k,j)
+ enddo
+ enddo
+ else
+#if (defined(AIX) && !defined(JUBL))
+ call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
+ if (iret.eq.0) goto 101
+ call xdrfint_(ixdrf, nss, iret)
+ if (iret.eq.0) goto 101
+ do j=1,nss
+ call xdrfint_(ixdrf, ihpb(j), iret)
+ if (iret.eq.0) goto 101
+ call xdrfint_(ixdrf, jhpb(j), iret)
+ if (iret.eq.0) goto 101
+ enddo
+ call xdrffloat_(ixdrf,reini,iret)
+ if (iret.eq.0) goto 101
+ call xdrffloat_(ixdrf,refree,iret)
+ if (iret.eq.0) goto 101
+ call xdrffloat_(ixdrf,rmsdev,iret)
+ if (iret.eq.0) goto 101
+ call xdrfint_(ixdrf,iscor,iret)
+ if (iret.eq.0) goto 101
+#else
+! write (iout,*) "calling xdrf3dfcoord"
+ call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
+! write (iout,*) "iret",iret
+! call flush(iout)
+ if (iret.eq.0) goto 101
+ call xdrfint(ixdrf, nss, iret)
+! write (iout,*) "iret",iret
+! write (iout,*) "nss",nss
+ call flush(iout)
+ if (iret.eq.0) goto 101
+ do k=1,nss
+ call xdrfint(ixdrf, ihpb(k), iret)
+ if (iret.eq.0) goto 101
+ call xdrfint(ixdrf, jhpb(k), iret)
+ if (iret.eq.0) goto 101
+ enddo
+ call xdrffloat(ixdrf,reini,iret)
+ if (iret.eq.0) goto 101
+ call xdrffloat(ixdrf,refree,iret)
+ if (iret.eq.0) goto 101
+ call xdrffloat(ixdrf,rmsdev,iret)
+ if (iret.eq.0) goto 101
+ call xdrfint(ixdrf,iscor,iret)
+ if (iret.eq.0) goto 101
+#endif
+ energy(jj+1)=reini
+ entfac(jj+1)=refree
+ rmstb(jj+1)=rmsdev
+ do k=1,nres
+ do l=1,3
+ c(l,k)=csingle(l,k)
+ enddo
+ enddo
+ do k=nnt,nct
+ do l=1,3
+ c(l,nres+k)=csingle(l,nres+k-nnt+1)
+ enddo
+ enddo
+ endif
+#ifdef DEBUG
+ write (iout,'(5hREAD ,i5,3f15.4,i10)') &
+ jj+1,energy(jj+1),entfac(jj+1),&
+ rmstb(jj+1),iscor
+ write (iout,*) "Conformation",jjj+1,jj+1
+ write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
+ write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
+ call flush(iout)
+#endif
+ call add_new_cconf(jjj,jj,jj_old,icount,Next)
+ enddo
+ 101 continue
+ write (iout,*) i-1," conformations read from DA file ",&
+ intinname(:ilen(intinname))
+ write (iout,*) jj," conformations read so far"
+ if (from_bx) then
+ close(intin)
+ else
+#if (defined(AIX) && !defined(JUBL))
+ call xdrfclose_(ixdrf, iret)
+#else
+ call xdrfclose(ixdrf, iret)
+#endif
+ endif
+#ifdef MPI
+!#ifdef DEBUG
+ write (iout,*) "jj_old",jj_old," jj",jj
+!#endif
+ call write_and_send_cconf(icount,jj_old,jj,Next)
+ call MPI_Send(0,1,MPI_INTEGER,Next,570,&
+ MPI_COMM_WORLD,IERROR)
+ jj_old=jj+1
+#else
+ call write_and_send_cconf(icount,jj_old,jj,Next)
+#endif
+ t_acq = tcpu() - t_acq
+#ifdef MPI
+ write (iout,*) "Processor",me,&
+ " time for conformation read/send",t_acq
+ ELSE
+! A worker gets the confs from the master and sends them to its neighbor
+ t_acq = tcpu()
+ call receive_and_pass_cconf(icount,jj_old,jj,&
+ Previous,Next)
+ t_acq = tcpu() - t_acq
+ ENDIF
+#endif
+ ncon=jj
+! close(icbase)
+ close(intin)
+
+ write(iout,*)"A total of",ncon," conformations read."
+
+ allocate(enetb(1:max_ene,ncon)) !(max_ene,maxstr_proc)
+#ifdef MPI
+! Check if everyone has the same number of conformations
+ call MPI_Allgather(ncon,1,MPI_INTEGER,&
+ ntot_all(0),1,MPI_INTEGER,MPI_Comm_World,IERROR)
+ lerr=.false.
+ do i=0,nprocs-1
+ if (i.ne.me) then
+ if (ncon.ne.ntot_all(i)) then
+ write (iout,*) "Number of conformations at processor",i,&
+ " differs from that at processor",me,&
+ ncon,ntot_all(i)
+ lerr = .true.
+ endif
+ endif
+ enddo
+ if (lerr) then
+ write (iout,*)
+ write (iout,*) "Number of conformations read by processors"
+ write (iout,*)
+ do i=0,nprocs-1
+ write (iout,'(8i10)') i,ntot_all(i)
+ enddo
+ write (iout,*) "Calculation terminated."
+ call flush(iout)
+ return 1
+ endif
+ return
+#endif
+ 1111 write(iout,*) "Error opening coordinate file ",&
+ intinname(:ilen(intinname))
+ call flush(iout)
+ return 1
+ end subroutine read_coords
+!------------------------------------------------------------------------------
+ subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
+
+ use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
+ use energy_data, only: itel,itype,dsc,max_ene
+ use control_data, only: symetr
+ use geometry, only: int_from_cart1
+! implicit none
+! include "DIMENSIONS"
+! include "sizesclu.dat"
+! include "COMMON.CLUSTER"
+! include "COMMON.CONTROL"
+! include "COMMON.CHAIN"
+! include "COMMON.INTERACT"
+! include "COMMON.LOCAL"
+! include "COMMON.IOUNITS"
+! include "COMMON.NAMES"
+! include "COMMON.VAR"
+! include "COMMON.SBRIDGE"
+! include "COMMON.GEO"
+ integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
+ nn,nn1,inan,Next,itj,chalen
+ real(kind=8) :: etot,energia(0:max_ene)
+ jjj=jjj+1
+ chalen=int((nct-nnt+2)/symetr)
+ call int_from_cart1(.false.)
+ do j=nnt+1,nct
+ if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
+ if (j.gt.2) then
+ if (itel(j).ne.0 .and. itel(j-1).ne.0) then
+ write (iout,*) "Conformation",jjj,jj+1
+ write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
+ chalen
+ write (iout,*) "The Cartesian geometry is:"
+ write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+ write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+ write (iout,*) "The internal geometry is:"
+ write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+ write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+ write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+ write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+ write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+ write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+ write (iout,*) &
+ "This conformation WILL NOT be added to the database."
+ return
+ endif
+ endif
+ endif
+ enddo
+ do j=nnt,nct
+ itj=itype(j)
+ if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))) &
+ .gt.2.0d0) then
+ write (iout,*) "Conformation",jjj,jj+1
+ write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
+ write (iout,*) "The Cartesian geometry is:"
+ write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+ write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+ write (iout,*) "The internal geometry is:"
+ write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+ write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+ write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+ write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+ write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+ write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+ write (iout,*) &
+ "This conformation WILL NOT be added to the database."
+ return
+ endif
+ enddo
+ do j=3,nres
+ if (theta(j).le.0.0d0) then
+ write (iout,*) &
+ "Zero theta angle(s) in conformation",jjj,jj+1
+ write (iout,*) "The Cartesian geometry is:"
+ write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
+ write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+ write (iout,*) "The internal geometry is:"
+ write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+ write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+ write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+ write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+ write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+ write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+ write (iout,*) &
+ "This conformation WILL NOT be added to the database."
+ return
+ endif
+ if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
+ enddo
+ jj=jj+1
+#ifdef DEBUG
+ write (iout,*) "Conformation",jjj,jj
+ write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
+ write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
+ write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+ write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+ write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+ write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+ write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
+ write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+ write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+ write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
+ write (iout,'(e15.5,16i5)') entfac(icount+1)
+! & iscore(icount+1,0)
+#endif
+ icount=icount+1
+ call store_cconf_from_file(jj,icount)
+ if (icount.eq.maxstr_proc) then
+#ifdef DEBUG
+ write (iout,* ) "jj_old",jj_old," jj",jj
+#endif
+ call write_and_send_cconf(icount,jj_old,jj,Next)
+ jj_old=jj+1
+ icount=0
+ endif
+ return
+ end subroutine add_new_cconf
+!------------------------------------------------------------------------------
+ subroutine store_cconf_from_file(jj,icount)
+
+ use energy_data, only: ihpb,jhpb
+! implicit none
+! include "DIMENSIONS"
+! include "sizesclu.dat"
+! include "COMMON.CLUSTER"
+! include "COMMON.CHAIN"
+! include "COMMON.SBRIDGE"
+! include "COMMON.INTERACT"
+! include "COMMON.IOUNITS"
+! include "COMMON.VAR"
+ integer :: i,j,jj,icount
+! Store the conformation that has been read in
+ do i=1,2*nres
+ do j=1,3
+ allcart(j,i,icount)=c(j,i)
+ enddo
+ enddo
+ nss_all(icount)=nss
+ do i=1,nss
+ ihpb_all(i,icount)=ihpb(i)
+ jhpb_all(i,icount)=jhpb(i)
+ enddo
+ return
+ end subroutine store_cconf_from_file
+!------------------------------------------------------------------------------
+ subroutine write_and_send_cconf(icount,jj_old,jj,Next)
+
+! implicit none
+! include "DIMENSIONS"
+! include "sizesclu.dat"
+#ifdef MPI
+ use MPI_data
+ include "mpif.h"
+ integer :: IERROR
+! include "COMMON.MPI"
+#endif
+! include "COMMON.CHAIN"
+! include "COMMON.SBRIDGE"
+! include "COMMON.INTERACT"
+! include "COMMON.IOUNITS"
+! include "COMMON.CLUSTER"
+! include "COMMON.VAR"
+ integer :: icount,jj_old,jj,Next
+! Write the structures to a scratch file
+#ifdef MPI
+! Master sends the portion of conformations that have been read in to the neighbor
+#ifdef DEBUG
+ write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
+ call flush(iout)
+#endif
+ call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
+ call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
+ Next,571,MPI_COMM_WORLD,IERROR)
+ call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
+ Next,572,MPI_COMM_WORLD,IERROR)
+ call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
+ Next,573,MPI_COMM_WORLD,IERROR)
+ call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
+ Next,577,MPI_COMM_WORLD,IERROR)
+ call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
+ Next,579,MPI_COMM_WORLD,IERROR)
+ call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
+ MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
+#endif
+ call dawrite_ccoords(jj_old,jj,icbase)
+ return
+ end subroutine write_and_send_cconf
+!------------------------------------------------------------------------------
+#ifdef MPI
+ subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
+
+ use MPI_data
+! implicit none
+! include "DIMENSIONS"
+! include "sizesclu.dat"
+ include "mpif.h"
+ integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
+! include "COMMON.MPI"
+! include "COMMON.CHAIN"
+! include "COMMON.SBRIDGE"
+! include "COMMON.INTERACT"
+! include "COMMON.IOUNITS"
+! include "COMMON.VAR"
+! include "COMMON.GEO"
+! include "COMMON.CLUSTER"
+ integer :: i,j,k,icount,jj_old,jj,Previous,Next
+ icount=1
+#ifdef DEBUG
+ write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
+ call flush(iout)
+#endif
+ do while (icount.gt.0)
+ call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
+ STATUS,IERROR)
+ call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
+ IERROR)
+#ifdef DEBUG
+ write (iout,*) "Processor",me," icount",icount
+#endif
+ if (icount.eq.0) return
+ call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
+ Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
+ call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
+ Next,571,MPI_COMM_WORLD,IERROR)
+ call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
+ Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
+ call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
+ Next,572,MPI_COMM_WORLD,IERROR)
+ call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
+ Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
+ call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
+ Next,573,MPI_COMM_WORLD,IERROR)
+ call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
+ Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
+ call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
+ Next,577,MPI_COMM_WORLD,IERROR)
+ call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
+ Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
+ call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
+ Next,579,MPI_COMM_WORLD,IERROR)
+ call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
+ MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
+ call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
+ MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
+ jj=jj_old+icount-1
+ call dawrite_ccoords(jj_old,jj,icbase)
+ jj_old=jj+1
+#ifdef DEBUG
+ write (iout,*) "Processor",me," received",icount," conformations"
+ do i=1,icount
+ write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
+ write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
+ write (iout,'(e15.5,16i5)') entfac(i)
+ enddo
+#endif
+ enddo
+ return
+ end subroutine receive_and_pass_cconf
+#endif
+!------------------------------------------------------------------------------
+ subroutine daread_ccoords(istart_conf,iend_conf)
+
+! implicit none
+! include "DIMENSIONS"
+! include "sizesclu.dat"
+#ifdef MPI
+ use MPI_data
+ include "mpif.h"
+! include "COMMON.MPI"
+#endif
+! include "COMMON.CHAIN"
+! include "COMMON.CLUSTER"
+! include "COMMON.IOUNITS"
+! include "COMMON.INTERACT"
+! include "COMMON.VAR"
+! include "COMMON.SBRIDGE"
+! include "COMMON.GEO"
+ integer :: istart_conf,iend_conf
+ integer :: i,j,ij,ii,iii
+ integer :: len
+ character(len=16) :: form,acc
+ character(len=32) :: nam
+!
+! Read conformations off a DA scratchfile.
+!
+#ifdef DEBUG
+ write (iout,*) "DAREAD_COORDS"
+ write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
+ inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
+ write (iout,*) "len=",len," form=",form," acc=",acc
+ write (iout,*) "nam=",nam
+ call flush(iout)
+#endif
+ do ii=istart_conf,iend_conf
+ ij = ii - istart_conf + 1
+ iii=list_conf(ii)
+#ifdef DEBUG
+ write (iout,*) "Reading binary file, record",iii," ii",ii
+ call flush(iout)
+#endif
+ read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
+ ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+ nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
+ entfac(ii),rmstb(ii)
+#ifdef DEBUG
+ write (iout,*) ii,iii,ij,entfac(ii)
+ write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
+ write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
+ i=nnt+nres,nct+nres)
+ write (iout,'(2e15.5)') entfac(ij)
+ write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
+ jhpb_all(i,ij),i=1,nss)
+ call flush(iout)
+#endif
+ enddo
+ return
+ end subroutine daread_ccoords
+!------------------------------------------------------------------------------
+ subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
+
+! implicit none
+! include "DIMENSIONS"
+! include "sizesclu.dat"
+#ifdef MPI
+ use MPI_data
+ include "mpif.h"
+! include "COMMON.MPI"
+#endif
+! include "COMMON.CHAIN"
+! include "COMMON.INTERACT"
+! include "COMMON.IOUNITS"
+! include "COMMON.VAR"
+! include "COMMON.SBRIDGE"
+! include "COMMON.GEO"
+! include "COMMON.CLUSTER"
+ integer :: istart_conf,iend_conf
+ integer :: i,j,ii,ij,iii,unit_out
+ integer :: len
+ character(len=16) :: form,acc
+ character(len=32) :: nam
+!
+! Write conformations to a DA scratchfile.
+!
+#ifdef DEBUG
+ write (iout,*) "DAWRITE_COORDS"
+ write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
+ write (iout,*) "lenrec",lenrec
+ inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
+ write (iout,*) "len=",len," form=",form," acc=",acc
+ write (iout,*) "nam=",nam
+ call flush(iout)
+#endif
+ do ii=istart_conf,iend_conf
+ iii=list_conf(ii)
+ ij = ii - istart_conf + 1
+#ifdef DEBUG
+ write (iout,*) "Writing binary file, record",iii," ii",ii
+ call flush(iout)
+#endif
+ write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
+ ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+ nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
+ entfac(ii),rmstb(ii)
+#ifdef DEBUG
+ write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
+ write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
+ nct+nres)
+ write (iout,'(2e15.5)') entfac(ij)
+ write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
+ nss_all(ij))
+ call flush(iout)
+#endif
+ enddo
+ return
+ end subroutine dawrite_ccoords
+!-----------------------------------------------------------------------------
+! readrtns.F
+!-----------------------------------------------------------------------------
+ subroutine read_control
+!
+! Read molecular data
+!
+ use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
+ use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
+ iscode,symetr,punch_dist,print_dist,nstart,nend,&
+ caonly,iopt,efree,lprint_cart,lprint_int
+! implicit none
+! include 'DIMENSIONS'
+! include 'sizesclu.dat'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CLUSTER'
+! include 'COMMON.CHAIN'
+! include 'COMMON.HEADER'
+! include 'COMMON.FFIELD'
+! include 'COMMON.FREE'
+! include 'COMMON.INTERACT'
+ character(len=320) :: controlcard !,ucase
+!#ifdef MPL
+! include 'COMMON.INFO'
+!#endif
+ integer :: i
+
+ read (INP,'(a80)') titel
+ call card_concat(controlcard,.true.)
+
+ call readi(controlcard,'NRES',nres,0)
+
+! call alloc_clust_arrays
+ allocate(rcutoff(max_cut+1)) !(max_cut+1)
+ allocate(beta_h(maxT)) !(maxT)
+ allocate(mult(nres)) !(maxres)
+
+
+ call readi(controlcard,'RESCALE',rescale_mode,2)
+ call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
+ write (iout,*) "DISTCHAINMAX",distchainmax
+ call readi(controlcard,'PDBOUT',outpdb,0)
+ call readi(controlcard,'MOL2OUT',outmol2,0)
+ refstr=(index(controlcard,'REFSTR').gt.0)
+ write (iout,*) "REFSTR",refstr
+ pdbref=(index(controlcard,'PDBREF').gt.0)
+ iscode=index(controlcard,'ONE_LETTER')
+ tree=(index(controlcard,'MAKE_TREE').gt.0)
+ min_var=(index(controlcard,'MINVAR').gt.0)
+ plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
+ punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
+ call readi(controlcard,'NCUT',ncut,1)
+ call readi(controlcard,'SYM',symetr,1)
+ write (iout,*) 'sym', symetr
+ call readi(controlcard,'NSTART',nstart,0)
+ call readi(controlcard,'NEND',nend,0)
+ call reada(controlcard,'ECUT',ecut,10.0d0)
+ call reada(controlcard,'PROB',prob_limit,0.99d0)
+ write (iout,*) "Probability limit",prob_limit
+ lgrp=(index(controlcard,'LGRP').gt.0)
+ caonly=(index(controlcard,'CA_ONLY').gt.0)
+ print_dist=(index(controlcard,'PRINT_DIST').gt.0)
+ call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
+ call readi(controlcard,'IOPT',iopt,2)
+ lside = index(controlcard,"SIDE").gt.0
+ efree = index(controlcard,"EFREE").gt.0
+ call readi(controlcard,'NTEMP',nT,1)
+ write (iout,*) "nT",nT
+!el call reada(controlcard,'TEMP0',temp0,300.0d0) !el
+ call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
+ write (iout,*) "nT",nT
+ write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
+ do i=1,nT
+ beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
+ enddo
+ write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
+ lprint_cart=index(controlcard,"PRINT_CART") .gt.0
+ lprint_int=index(controlcard,"PRINT_INT") .gt.0
+ if (min_var) iopt=1
+ return
+ end subroutine read_control
+!-----------------------------------------------------------------------------
+ subroutine molread
+!
+! Read molecular data.
+!
+ use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
+ use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
+! wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
+! wturn3,wturn4,wturn6,wvdwpp,weights
+ use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
+ indpdb
+ use geometry, only: chainbuild,alloc_geo_arrays
+ use energy, only: alloc_ener_arrays
+ use control, only: rescode,setup_var,init_int_table
+ use conform_compar, only: contact
+! implicit none
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.LOCAL'
+! include 'COMMON.NAMES'
+! include 'COMMON.CHAIN'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.HEADER'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CONTACTS'
+! include 'COMMON.TIME1'
+!#ifdef MPL
+! include 'COMMON.INFO'
+!#endif
+ character(len=4) :: sequence(nres) !(maxres)
+ character(len=800) :: weightcard
+! integer rescode
+ real(kind=8) :: x(6*nres) !(maxvar)
+ integer :: itype_pdb(nres) !(maxres)
+! logical seq_comp
+ integer :: i,j,kkk
+!
+! Body
+!
+!el allocate(weights(n_ene))
+ allocate(weights(max_ene))
+ call alloc_geo_arrays
+ call alloc_ener_arrays
+!-----------------------------
+ allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
+ allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
+ allocate(itype(nres+2)) !(maxres)
+ allocate(itel(nres+2))
+
+ do i=1,2*nres+2
+ do j=1,3
+ c(j,i)=0
+ dc(j,i)=0
+ enddo
+ enddo
+ do i=1,nres+2
+ itype(i)=0
+ itel(i)=0
+ enddo
+!--------------------------
+! Read weights of the subsequent energy terms.
+ call card_concat(weightcard,.true.)
+ call reada(weightcard,'WSC',wsc,1.0d0)
+ call reada(weightcard,'WLONG',wsc,wsc)
+ call reada(weightcard,'WSCP',wscp,1.0d0)
+ call reada(weightcard,'WELEC',welec,1.0D0)
+ call reada(weightcard,'WVDWPP',wvdwpp,welec)
+ call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
+ call reada(weightcard,'WCORR4',wcorr4,0.0D0)
+ call reada(weightcard,'WCORR5',wcorr5,0.0D0)
+ call reada(weightcard,'WCORR6',wcorr6,0.0D0)
+ call reada(weightcard,'WTURN3',wturn3,1.0D0)
+ call reada(weightcard,'WTURN4',wturn4,1.0D0)
+ call reada(weightcard,'WTURN6',wturn6,1.0D0)
+ call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
+ call reada(weightcard,'WSCCOR',wsccor,1.0D0)
+ call reada(weightcard,'WBOND',wbond,1.0D0)
+ call reada(weightcard,'WTOR',wtor,1.0D0)
+ call reada(weightcard,'WTORD',wtor_d,1.0D0)
+ call reada(weightcard,'WANG',wang,1.0D0)
+ call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+ call reada(weightcard,'SCAL14',scal14,0.4D0)
+ call reada(weightcard,'SCALSCP',scalscp,1.0d0)
+ call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
+ call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
+ call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el
+ if (index(weightcard,'SOFT').gt.0) ipot=6
+! 12/1/95 Added weight for the multi-body term WCORR
+ call reada(weightcard,'WCORRH',wcorr,1.0D0)
+ if (wcorr4.gt.0.0d0) wcorr=wcorr4
+ weights(1)=wsc
+ weights(2)=wscp
+ weights(3)=welec
+ weights(4)=wcorr
+ weights(5)=wcorr5
+ weights(6)=wcorr6
+ weights(7)=wel_loc
+ weights(8)=wturn3
+ weights(9)=wturn4
+ weights(10)=wturn6
+ weights(11)=wang
+ weights(12)=wscloc
+ weights(13)=wtor
+ weights(14)=wtor_d
+ weights(15)=wstrain
+ weights(16)=wvdwpp
+ weights(17)=wbond
+ weights(18)=scal14
+!el weights(19)=wsccor !!!!!!!!!!!!!!!!
+ weights(21)=wsccor
+ write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
+ wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
+ wturn4,wturn6,wsccor
+ 10 format (/'Energy-term weights (unscaled):'// &
+ 'WSCC= ',f10.6,' (SC-SC)'/ &
+ 'WSCP= ',f10.6,' (SC-p)'/ &
+ 'WELEC= ',f10.6,' (p-p electr)'/ &
+ 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
+ 'WBOND= ',f10.6,' (stretching)'/ &
+ 'WANG= ',f10.6,' (bending)'/ &
+ 'WSCLOC= ',f10.6,' (SC local)'/ &
+ 'WTOR= ',f10.6,' (torsional)'/ &
+ 'WTORD= ',f10.6,' (double torsional)'/ &
+ 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
+ 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
+ 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
+ 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
+ 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
+ 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
+ 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
+ 'WTURN6= ',f10.6,' (turns, 6th order)'/ &
+ 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
+
+ if (wcorr4.gt.0.0d0) then
+ write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
+ 'between contact pairs of peptide groups'
+ write (iout,'(2(a,f5.3/))') &
+ 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
+ 'Range of quenching the correlation terms:',2*delt_corr
+ else if (wcorr.gt.0.0d0) then
+ write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
+ 'between contact pairs of peptide groups'
+ endif
+ write (iout,'(a,f8.3)') &
+ 'Scaling factor of 1,4 SC-p interactions:',scal14
+ write (iout,'(a,f8.3)') &
+ 'General scaling factor of SC-p interactions:',scalscp
+ r0_corr=cutoff_corr-delt_corr
+ do i=1,20
+ aad(i,1)=scalscp*aad(i,1)
+ aad(i,2)=scalscp*aad(i,2)
+ bad(i,1)=scalscp*bad(i,1)
+ bad(i,2)=scalscp*bad(i,2)
+ enddo
+
+ call flush(iout)
+ print *,'indpdb=',indpdb,' pdbref=',pdbref
+
+! Read sequence if not taken from the pdb file.
+ if (iscode.gt.0) then
+ read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
+ else
+ read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
+ endif
+! Convert sequence to numeric code
+ do i=1,nres
+ itype(i)=rescode(i,sequence(i),iscode)
+ enddo
+ print *,nres
+ print '(20i4)',(itype(i),i=1,nres)
+
+ do i=1,nres
+#ifdef PROCOR
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
+#else
+ if (itype(i).eq.ntyp1) then
+#endif
+ itel(i)=0
+#ifdef PROCOR
+ else if (iabs(itype(i+1)).ne.20) then
+#else
+ else if (iabs(itype(i)).ne.20) then
+#endif
+ itel(i)=1
+ else
+ itel(i)=2
+ endif
+ enddo
+ write (iout,*) "ITEL"
+ do i=1,nres-1
+ write (iout,*) i,itype(i),itel(i)
+ enddo
+
+ print *,'Call Read_Bridge.'
+ call read_bridge
+ nnt=1
+ nct=nres
+ print *,'NNT=',NNT,' NCT=',NCT
+ if (itype(1).eq.ntyp1) nnt=2
+ if (itype(nres).eq.ntyp1) nct=nct-1
+ if (nstart.lt.nnt) nstart=nnt
+ if (nend.gt.nct .or. nend.eq.0) nend=nct
+ write (iout,*) "nstart",nstart," nend",nend
+ nres0=nres
+! if (pdbref) then
+! read(inp,'(a)') pdbfile
+! write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
+! open(ipdbin,file=pdbfile,status='old',err=33)
+! goto 34
+! 33 write (iout,'(a)') 'Error opening PDB file.'
+! stop
+! 34 continue
+! print *,'Begin reading pdb data'
+! call readpdb
+! print *,'Finished reading pdb data'
+! write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
+! do i=1,nres
+! itype_pdb(i)=itype(i)
+! enddo
+! close (ipdbin)
+! write (iout,'(a,i3)') 'nsup=',nsup
+! nstart_seq=nnt
+! if (nsup.le.(nct-nnt+1)) then
+! do i=0,nct-nnt+1-nsup
+! if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
+! nstart_seq=nnt+i
+! goto 111
+! endif
+! enddo
+! write (iout,'(a)')
+! & 'Error - sequences to be superposed do not match.'
+! stop
+! else
+! do i=0,nsup-(nct-nnt+1)
+! if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
+! & then
+! nstart_sup=nstart_sup+i
+! nsup=nct-nnt+1
+! goto 111
+! endif
+! enddo
+! write (iout,'(a)')
+! & 'Error - sequences to be superposed do not match.'
+! endif
+! 111 continue
+! write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
+! & ' nstart_seq=',nstart_seq
+! endif
+write(iout,*)"przed ini_int_tab"
+ call init_int_table
+write(iout,*)"po ini_int_tab"
+write(iout,*)"przed setup var"
+ call setup_var
+write(iout,*)"po setup var"
+ write (iout,*) "molread: REFSTR",refstr
+ if (refstr) then
+ if (.not.pdbref) then
+ call read_angles(inp,*38)
+ goto 39
+ 38 write (iout,'(a)') 'Error reading reference structure.'
+#ifdef MPL
+ call mp_stopall(Error_Msg)
+#else
+ stop 'Error reading reference structure'
+#endif
+ 39 call chainbuild
+ nstart_sup=nnt
+ nstart_seq=nnt
+ nsup=nct-nnt+1
+ kkk=1
+ do i=1,2*nres
+ do j=1,3
+ cref(j,i,kkk)=c(j,i)
+ enddo
+ enddo
+ endif
+ call contact(.true.,ncont_ref,icont_ref)
+ endif
+ return
+ end subroutine molread
+!-----------------------------------------------------------------------------
+ subroutine openunits
+! implicit none
+! include 'DIMENSIONS'
+ use control_data, only: from_cx,from_bx,from_cart
+#ifdef MPI
+ use MPI_data
+ include "mpif.h"
+ character(len=3) :: liczba
+! include "COMMON.MPI"
+#endif
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CONTROL'
+ integer :: lenpre,lenpot !,ilen
+! external ilen
+ character(len=16) :: cformat,cprint
+! character(len=16) ucase
+ integer :: lenint,lenout
+ call getenv('INPUT',prefix)
+ call getenv('OUTPUT',prefout)
+ call getenv('INTIN',prefintin)
+ call getenv('COORD',cformat)
+ call getenv('PRINTCOOR',cprint)
+ call getenv('SCRATCHDIR',scratchdir)
+ from_bx=.true.
+ from_cx=.false.
+ if (index(ucase(cformat),'CX').gt.0) then
+ from_cx=.true.
+ from_bx=.false.
+ endif
+ from_cart=.true.
+ lenpre=ilen(prefix)
+ lenout=ilen(prefout)
+ lenint=ilen(prefintin)
+! Get the names and open the input files
+ open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
+#ifdef MPI
+ write (liczba,'(bz,i3.3)') me
+ outname=prefout(:lenout)//'_clust.out_'//liczba
+#else
+ outname=prefout(:lenout)//'_clust.out'
+#endif
+ if (from_bx) then
+ intinname=prefintin(:lenint)//'.bx'
+ else if (from_cx) then
+ intinname=prefintin(:lenint)//'.cx'
+ else
+ intinname=prefintin(:lenint)//'.int'
+ endif
+ rmsname=prefintin(:lenint)//'.rms'
+ open (jplot,file=prefout(:ilen(prefout))//'.tex',&
+ status='unknown')
+ open (jrms,file=rmsname,status='unknown')
+ open(iout,file=outname,status='unknown')
+! Get parameter filenames and open the parameter files.
+ call getenv('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old')
+ call getenv('THETPAR',thetname)
+ open (ithep,file=thetname,status='old')
+ call getenv('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old')
+ call getenv('TORPAR',torname)
+ open (itorp,file=torname,status='old')
+ call getenv('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old')
+ call getenv('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old')
+ call getenv('ELEPAR',elename)
+ open (ielep,file=elename,status='old')
+ call getenv('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old')
+ call getenv('SIDEP',sidepname)
+ open (isidep1,file=sidepname,status="old")
+ call getenv('SCCORPAR',sccorname)
+ open (isccor,file=sccorname,status="old")
+#ifndef OLDSCP
+!
+! 8/9/01 In the newest version SCp interaction constants are read from a file
+! Use -DOLDSCP to use hard-coded constants instead.
+!
+ call getenv('SCPPAR',scpname)
+ open (iscpp,file=scpname,status='old')
+#endif
+ return
+ end subroutine openunits
+!-----------------------------------------------------------------------------
+! geomout.F
+!-----------------------------------------------------------------------------
+ subroutine pdboutC(etot,rmsd,tytul)
+
+ use energy_data, only: ihpb,jhpb,itype
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.HEADER'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.TEMPFAC'
+ character(len=50) :: tytul
+ character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
+ 'G','H','I','J'/)
+ integer :: ica(nres)
+ real(kind=8) :: etot,rmsd
+ integer :: iatom,ichain,ires,i,j,iti
+
+ write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
+ ' ENERGY ',etot,' RMS ',rmsd
+ iatom=0
+ ichain=1
+ ires=0
+ do i=nnt,nct
+ iti=itype(i)
+ if (iti.eq.ntyp1) then
+ ichain=ichain+1
+ ires=0
+ write (ipdb,'(a)') 'TER'
+ else
+ ires=ires+1
+ iatom=iatom+1
+ ica(i)=iatom
+ write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
+ ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
+ if (iti.ne.10) then
+ iatom=iatom+1
+ write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
+ ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
+ endif
+ endif
+ enddo
+ write (ipdb,'(a)') 'TER'
+ do i=nnt,nct-1
+ if (itype(i).eq.ntyp1) cycle
+ if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
+ write (ipdb,30) ica(i),ica(i+1)
+ else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+ write (ipdb,30) ica(i),ica(i+1),ica(i)+1
+ else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
+ write (ipdb,30) ica(i),ica(i)+1
+ endif
+ enddo
+ if (itype(nct).ne.10) then
+ write (ipdb,30) ica(nct),ica(nct)+1
+ endif
+ do i=1,nss
+ write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
+ enddo
+ write (ipdb,'(a6)') 'ENDMDL'
+ 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
+ 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
+ 30 FORMAT ('CONECT',8I5)
+ return
+ end subroutine pdboutC
+!-----------------------------------------------------------------------------
+ subroutine cartout(igr,i,etot,free,rmsd,plik)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'sizesclu.dat'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.GEO'
+! include 'COMMON.CLUSTER'
+ integer :: igr,i,j,k
+ real(kind=8) :: etot,free,rmsd
+ character(len=80) :: plik
+ open (igeom,file=plik,position='append')
+ write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
+ write (igeom,'(i4,$)') &
+ nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
+ write (igeom,'(i10)') iscore(i)
+ write (igeom,'(8f10.5)') &
+ ((allcart(k,j,i),k=1,3),j=1,nres),&
+ ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
+ return
+ end subroutine cartout
+!------------------------------------------------------------------------------
+! subroutine alloc_clust_arrays(n_conf)
+
+! integer :: n_conf
+!COMMON.CLUSTER
+! common /clu/
+! allocate(diss(maxdist)) !(maxdist)
+!el allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
+! allocatable :: enetb !(max_ene,maxstr_proc)
+!el allocate(entfac(maxconf)) !(maxconf)
+! allocatable :: totfree_gr !(maxgr)
+!el allocate(rcutoff(max_cut+1)) !(max_cut+1)
+! common /clu1/
+! allocatable :: licz,iass !(maxgr)
+! allocatable :: nconf !(maxgr,maxingr)
+! allocatable :: iass_tot !(maxgr,max_cut)
+! allocatable :: list_conf !(maxconf)
+! common /alles/
+!el allocatable :: allcart !(3,maxres2,maxstr_proc)
+!el allocate(rmstb(maxconf)) !(maxconf)
+!el allocate(mult(nres)) !(maxres)
+!el allocatable :: nss_all !(maxstr_proc)
+!el allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
+! allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
+!COMMON.TEMPFAC
+! common /factemp/
+! allocatable :: tempfac !(2,maxres)
+!COMMON.FREE
+! common /free/
+!el allocate(beta_h(maxT)) !(maxT)
+
+! end subroutine alloc_clust_arrays
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ end module io_clust
projects / unres4.git / blobdiff