module io_clust !----------------------------------------------------------------------------- use clust_data use io_units ! use names use io_base !, only: ilen use geometry_data, only: nres,c,nres_molec,boxxsize,boxysize,boxzsize 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,molnum 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,mnum 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 mnum=molnum(j) if (mnum.eq.5) cycle if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) 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 mnum=molnum(j) itj=itype(j,mnum) if (mnum.eq.5) cycle if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) & .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,rlamb_ele,& r_cut_ele ! 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_molec(1),0) call readi(controlcard,"NRES_NUCL",nres_molec(2),0) call readi(controlcard,"NRES_CAT",nres_molec(5),0) nres=0 do i=1,5 nres=nres_molec(i)+nres enddo ! 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 reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) ! ions=index(controlcard,"IONS").gt.0 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0) call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0) write(iout,*) "R_CUT_ELE=",r_cut_ele 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 io_base, only: rescode use control, only: 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,mnum ! ! 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,5)) !(maxres) allocate(molnum(nres+1)) 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 j=1,5 do i=1,nres+2 itype(i,j)=0 itel(i)=0 enddo 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,'WCATCAT',wcatcat,0.0d0) call reada(weightcard,'WCATPROT',wcatprot,0.0d0) 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_molec(1) mnum=1 molnum(i)=1 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum) enddo do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2) mnum=2 molnum(i)=2 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum) enddo do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5) mnum=5 molnum(i)=5 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum) enddo print *,nres print '(20i4)',(itype(i,molnum(i)),i=1,nres) do i=1,nres-1 #ifdef PROCOR if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))& .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then #else if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then #endif itel(i)=0 #ifdef PROCOR else if (iabs(itype(i+1,1)).ne.20) then #else else if (iabs(itype(i,1)).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,molnum(i)),itel(i) enddo print *,'Call Read_Bridge.' call read_bridge nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) 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") call getenv('THETPAR_NUCL',thetname_nucl) open (ithep_nucl,file=thetname_nucl,status='old') call getenv('ROTPAR_NUCL',rotname_nucl) open (irotam_nucl,file=rotname_nucl,status='old') call getenv('TORPAR_NUCL',torname_nucl) open (itorp_nucl,file=torname_nucl,status='old') call getenv('TORDPAR_NUCL',tordname_nucl) open (itordp_nucl,file=tordname_nucl,status='old') call getenv('SIDEPAR_NUCL',sidename_nucl) open (isidep_nucl,file=sidename_nucl,status='old') call getenv('SIDEPAR_SCBASE',sidename_scbase) open (isidep_scbase,file=sidename_scbase,status='old') call getenv('PEPPAR_PEPBASE',pepname_pepbase) open (isidep_pepbase,file=pepname_pepbase,status='old') call getenv('SCPAR_PHOSPH',pepname_scpho) open (isidep_scpho,file=pepname_scpho,status='old') call getenv('PEPPAR_PHOSPH',pepname_peppho) open (isidep_peppho,file=pepname_peppho,status='old') call getenv('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old') call getenv('TUBEPAR',tubename) open (itube,file=tubename,status='old') call getenv('IONPAR',ionname) open (iion,file=ionname,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,molnum ! 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,mnum,mnum1,boxxxshift(3) real(kind=8) :: boxxxx(3) write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),& ' ENERGY ',etot,' RMS ',rmsd iatom=0 ichain=1 ires=0 boxxxshift(1)=int(c(1,nnt)/boxxsize) ! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1 boxxxshift(2)=int(c(2,nnt)/boxzsize) ! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1 boxxxshift(3)=int(c(3,nnt)/boxzsize) ! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1 boxxxx(1)=boxxsize boxxxx(2)=boxysize boxxxx(3)=boxzsize do i=nnt,nct mnum=molnum(i) iti=itype(i,mnum) if (iti.eq.ntyp1_molec(mnum)) then ichain=ichain+1 ires=0 write (ipdb,'(a)') 'TER' else ires=ires+1 iatom=iatom+1 ica(i)=iatom if (mnum.eq.5) then do j=1,3 if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j) else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j) else c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j) endif enddo endif write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i) if ((iti.ne.10).and.(mnum.ne.5)) then iatom=iatom+1 write (ipdb,20) iatom,restyp(iti,mnum),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 mnum=molnum(i) mnum1=molnum(i+1) if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then write (ipdb,30) ica(i),ica(i+1) else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then write (ipdb,30) ica(i),ica(i+1),ica(i)+1 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then write (ipdb,30) ica(i),ica(i)+1 endif enddo if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) 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