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) nss_all=0 ihpb_all=0 jhpb_all=0 #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,i,"TUTU" 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 ! write(iout,*), ihpb(k),jhpb(k),"TUTU" enddo call xdrffloat(ixdrf,reini,iret) if (iret.eq.0) goto 101 ! write(iout,*) "TUTU", reini call xdrffloat(ixdrf,refree,iret) ! write(iout,*) "TUTU", refree 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-1 mnum=molnum(j) ! write (iout,*) "Check atom",j if (mnum.ne.1) 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.ne.1) 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) use energy_data, only: dyn_ss ! 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 if (dyn_ss) then 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),& entfac(ii),rmstb(ii) else 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) endif #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" use energy_data, only: dyn_ss #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 if (dyn_ss) then 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),& entfac(ii),rmstb(ii) else 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) endif #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,nclust,tor_mode,scelemode ! 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 print *,"TU",nres_molec(:) ! 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) call readi(controlcard,'NCLUST',nclust,5) ! ions=index(controlcard,"IONS").gt.0 call readi(controlcard,"SCELEMODE",scelemode,0) print *,"SCELE",scelemode call readi(controlcard,'TORMODE',tor_mode,0) !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "torsional and valence angle mode",tor_mode call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.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) real(kind=8) :: ssscale 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 call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0) call reada(weightcard,'WELPP',welpp,0.0d0) call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0) call reada(weightcard,'WELPSB',welpsb,0.0D0) call reada(weightcard,'WVDWSB',wvdwsb,0.0d0) call reada(weightcard,'WELSB',welsb,0.0D0) call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0) call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0) call reada(weightcard,'WSBLOC',wsbloc,0.0D0) call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0) ! print *,"KUR..",wtor_nucl call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0) call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0) call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0) ! 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) call reada(weightcard,'WSCBASE',wscbase,0.0D0) call reada(weightcard,'WPEPBASE',wpepbase,1.0d0) call reada(weightcard,'WSCPHO',wscpho,0.0d0) call reada(weightcard,'WPEPPHO',wpeppho,0.0d0) call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0) call reada(weightcard,"D0CM",d0cm,3.78d0) call reada(weightcard,"AKCM",akcm,15.1d0) call reada(weightcard,"AKTH",akth,11.0d0) call reada(weightcard,"AKCT",akct,12.0d0) call reada(weightcard,"V1SS",v1ss,-1.08d0) call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) call reada(weightcard,"ATRISS",atriss,0.301D0) call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) call reada(weightcard,"DTRISS",dtriss,1.001D0) call reada(weightcard,"SSSCALE",ssscale,1.0D0) dyn_ss=(index(weightcard,'DYN_SS').gt.0) call reada(weightcard,"HT",Ht,0.0D0) if (dyn_ss) then ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale Ht=(Ht/wsc-0.25*eps(1,1))*ssscale akcm=akcm/wsc*ssscale akth=akth/wsc*ssscale akct=akct/wsc*ssscale v1ss=v1ss/wsc*ssscale v2ss=v2ss/wsc*ssscale v3ss=v3ss/wsc*ssscale else ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain endif 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 weights(26)=wvdwpp_nucl weights(27)=welpp weights(28)=wvdwpsb weights(29)=welpsb weights(30)=wvdwsb weights(31)=welsb weights(32)=wbond_nucl weights(33)=wang_nucl weights(34)=wsbloc weights(35)=wtor_nucl weights(36)=wtor_d_nucl weights(37)=wcorr_nucl weights(38)=wcorr3_nucl weights(41)=wcatcat weights(42)=wcatprot weights(46)=wscbase weights(47)=wscpho weights(48)=wpeppho weights(49)=wpeppho weights(50)=wcatnucl 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 print *,nres_molec(:),nres ! 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 write (iout,*),i,sequence(i) 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" if (ns.gt.0.and.dyn_ss) then do i=nss+1,nhpb ihpb(i-nss)=ihpb(i) jhpb(i-nss)=jhpb(i) forcon(i-nss)=forcon(i) dhpb(i-nss)=dhpb(i) enddo nhpb=nhpb-nss nss=0 link_start=1 link_end=nhpb ! call hpb_partition do i=1,ns dyn_ss_mask(iss(i))=.true. enddo endif 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('BONDPAR_NUCL',bondname_nucl) open (ibond_nucl,file=bondname_nucl,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') call getenv('IONPAR_NUCL',ionnuclname) open (iionnucl,file=ionnuclname,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') call getenv('SCPPAR_NUCL',scpname_nucl) open (iscpp_nucl,file=scpname_nucl,status='old') #endif return end subroutine openunits !----------------------------------------------------------------------------- ! geomout.F !----------------------------------------------------------------------------- subroutine pdboutC(etot,rmsd,tytul) use energy_data, only: ihpb,jhpb,itype,molnum use energy, only:boxshift,to_box ! 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),k,iti1 real(kind=8) :: etot,rmsd integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3) real(kind=8) :: boxxxx(3),cbeg(3),Rdist(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 ! do i=1,3 ! if (boxxxshift(i).gt.2) boxxxshift(i)=2 ! if (boxxxshift(i).lt.-2) boxxxshift(i)=-2 ! ! enddo do j=1,3 cbeg(j)=c(j,nnt) enddo call to_box(cbeg(1),cbeg(2),cbeg(3)) boxxxx(1)=boxxsize boxxxx(2)=boxysize boxxxx(3)=boxzsize do i=nnt,nct mnum=molnum(i) iti=itype(i,mnum) if (i.ne.nct) then iti1=itype(i+1,mnum) mnum1=molnum(i+1) if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle endif if (iti.eq.ntyp1_molec(mnum)) then ichain=ichain+1 ires=0 write (ipdb,'(a)') 'TER' else ires=ires+1 if ((ires.eq.2).and.(mnum.ne.5)) then do j=1,3 cbeg(j)=c(j,i-1) enddo endif iatom=iatom+1 ica(i)=iatom call to_box(c(1,i),c(2,i),c(3,i)) DO k=1,3 Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k)) c(k,i)=cbeg(k)+Rdist(k) ENDDO if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres)) DO k=1,3 Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k)) c(k,i+nres)=cbeg(k)+Rdist(k) 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