X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fio_clust.f90;fp=source%2Fcluster%2Fio_clust.f90;h=db74bf627a814128bcef109a907719154df6f51e;hb=299e2c41124d3fa8adba7244716515a2cc160ed1;hp=0000000000000000000000000000000000000000;hpb=65b670d78bae7c2df18f3465fdad7b4f13490e33;p=unres4.git diff --git a/source/cluster/io_clust.f90 b/source/cluster/io_clust.f90 new file mode 100644 index 0000000..db74bf6 --- /dev/null +++ b/source/cluster/io_clust.f90 @@ -0,0 +1,1824 @@ + 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