module io !----------------------------------------------------------------------- use io_units use names use io_base use io_config implicit none !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !----------------------------------------------------------------------------- ! bank.F io_csa !----------------------------------------------------------------------------- subroutine write_csa_pdb(var,ene,nft,ik,iw_pdb) use csa_data use geometry_data, only:nres,nvar use geometry, only:var_to_geom,chainbuild use compare, only:secondary2 ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.VAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.MINIM' ! include 'COMMON.SETUP' ! include 'COMMON.GEO' ! include 'COMMON.CHAIN' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.SBRIDGE' integer :: lenpre,lenpot !,ilen !el external ilen real(kind=8),dimension(nvar) :: var !(maxvar) (maxvar=6*maxres) character(len=50) :: titelloc character(len=3) :: zahl real(kind=8),dimension(mxch*(mxch+1)/2+1) :: ene !el local variables integer :: nft,ik,iw_pdb nmin_csa=nmin_csa+1 if(ene(1).lt.eglob_csa) then eglob_csa=ene(1) nglob_csa=nglob_csa+1 call numstr(nglob_csa,zahl) call var_to_geom(nvar,var) call chainbuild call secondary2(.false.) lenpre=ilen(prefix) open(icsa_pdb,file=prefix(:lenpre)//'@'//zahl//'.pdb') if (iw_pdb.eq.1) then write(titelloc,'(a2,i3,a3,i9,a3,i6)') & 'GM',nglob_csa,' e ',nft,' m ',nmin_csa else write(titelloc,'(a2,i3,a3,i9,a3,i6,a5,f5.2,a5,f5.1)') & 'GM',nglob_csa,' e ',nft,' m ',nmin_csa,' rms ',& rmsn(ik),' %NC ',pncn(ik)*100 endif call pdbout(eglob_csa,titelloc,icsa_pdb) close(icsa_pdb) endif return end subroutine write_csa_pdb !----------------------------------------------------------------------------- ! csa.f io_csa !----------------------------------------------------------------------------- subroutine from_pdb(n,idum) ! This subroutine stores the UNRES int variables generated from ! subroutine readpdb into the 1st conformation of in dihang_in. ! Subsequent n-1 conformations of dihang_in have identical values ! of theta and phi as the 1st conformation but random values for ! alph and omeg. ! The array cref (also generated from subroutine readpdb) is stored ! to crefjlee to be used for rmsd calculation in CSA, if necessary. use csa_data use geometry_data use random, only: ran1 ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.BANK' ! include 'COMMON.GEO' !el local variables integer :: n,idum,m,i,j,k,kk,kkk real(kind=8) :: e m=1 do j=2,nres-1 dihang_in(1,j,1,m)=theta(j+1) dihang_in(2,j,1,m)=phi(j+2) dihang_in(3,j,1,m)=alph(j) dihang_in(4,j,1,m)=omeg(j) enddo dihang_in(2,nres-1,1,k)=0.0d0 do m=2,n do k=2,nres-1 dihang_in(1,k,1,m)=dihang_in(1,k,1,1) dihang_in(2,k,1,m)=dihang_in(2,k,1,1) if(dabs(dihang_in(3,k,1,1)).gt.1.d-6) then dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0 dihang_in(3,k,1,m)=dihang_in(3,k,1,m)*deg2rad endif if(dabs(dihang_in(4,k,1,1)).gt.1.d-6) then dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0 dihang_in(4,k,1,m)=dihang_in(4,k,1,m)*deg2rad endif enddo enddo ! Store cref to crefjlee (they are in COMMON.CHAIN). do k=1,2*nres do kk=1,3 kkk=1 crefjlee(kk,k)=cref(kk,k,kkk) enddo enddo open(icsa_native_int,file=csa_native_int,status="old") do m=1,n write(icsa_native_int,*) m,e write(icsa_native_int,200) & (dihang_in(1,k,1,m)*rad2deg,k=2,nres-1) write(icsa_native_int,200) & (dihang_in(2,k,1,m)*rad2deg,k=2,nres-2) write(icsa_native_int,200) & (dihang_in(3,k,1,m)*rad2deg,k=2,nres-1) write(icsa_native_int,200) & (dihang_in(4,k,1,m)*rad2deg,k=2,nres-1) enddo do k=1,nres write(icsa_native_int,200) (crefjlee(i,k),i=1,3) enddo close(icsa_native_int) 200 format (8f10.4) return end subroutine from_pdb !----------------------------------------------------------------------------- subroutine from_int(n,mm,idum) use csa_data use geometry_data use energy_data use geometry, only:chainbuild,gen_side use energy, only:etotal use compare ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.INTERACT' ! include 'COMMON.BANK' ! include 'COMMON.GEO' ! include 'COMMON.CONTACTS' ! integer ilen !el external ilen logical :: fail real(kind=8),dimension(0:n_ene) :: energia !el local variables integer :: n,mm,idum,i,ii,j,m,k,kk,maxcount_fail,icount_fail,maxsi real(kind=8) :: co open(icsa_native_int,file=csa_native_int,status="old") read (icsa_native_int,*) call read_angles(icsa_native_int,*10) goto 11 10 write (iout,'(2a)') "CHUJ NASTAPIL - error in ",& csa_native_int(:ilen(csa_native_int)) 11 continue call intout do j=2,nres-1 dihang_in(1,j,1,1)=theta(j+1) dihang_in(2,j,1,1)=phi(j+2) dihang_in(3,j,1,1)=alph(j) dihang_in(4,j,1,1)=omeg(j) enddo dihang_in(2,nres-1,1,1)=0.0d0 ! read(icsa_native_int,*) ind,e ! read(icsa_native_int,200) (dihang_in(1,k,1,1),k=2,nres-1) ! read(icsa_native_int,200) (dihang_in(2,k,1,1),k=2,nres-2) ! read(icsa_native_int,200) (dihang_in(3,k,1,1),k=2,nres-1) ! read(icsa_native_int,200) (dihang_in(4,k,1,1),k=2,nres-1) ! dihang_in(2,nres-1,1,1)=0.d0 maxsi=100 maxcount_fail=100 do m=mm+2,n ! do k=2,nres-1 ! dihang_in(1,k,1,m)=dihang_in(1,k,1,1) ! dihang_in(2,k,1,m)=dihang_in(2,k,1,1) ! if(abs(dihang_in(3,k,1,1)).gt.1.d-3) then ! dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0 ! endif ! if(abs(dihang_in(4,k,1,1)).gt.1.d-3) then ! dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0 ! endif ! enddo ! call intout fail=.true. icount_fail=0 DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL) do i=nnt,nct if (itype(i,1).ne.10) then !d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1) fail=.true. ii=0 do while (fail .and. ii .le. maxsi) call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,molnum(i)) ii = ii+1 enddo endif enddo call chainbuild call etotal(energia) fail = (energia(0).ge.1.0d20) icount_fail=icount_fail+1 ENDDO if (icount_fail.gt.maxcount_fail) then write (iout,*) & 'Failed to generate non-overlaping near-native conf.',& m endif do j=2,nres-1 dihang_in(1,j,1,m)=theta(j+1) dihang_in(2,j,1,m)=phi(j+2) dihang_in(3,j,1,m)=alph(j) dihang_in(4,j,1,m)=omeg(j) enddo dihang_in(2,nres-1,1,m)=0.0d0 enddo ! do m=1,n ! write(icsa_native_int,*) m,e ! write(icsa_native_int,200) (dihang_in(1,k,1,m),k=2,nres-1) ! write(icsa_native_int,200) (dihang_in(2,k,1,m),k=2,nres-2) ! write(icsa_native_int,200) (dihang_in(3,k,1,m),k=2,nres-1) ! write(icsa_native_int,200) (dihang_in(4,k,1,m),k=2,nres-1) ! enddo ! close(icsa_native_int) ! do m=mm+2,n ! do i=1,4 ! do j=2,nres-1 ! dihang_in(i,j,1,m)=dihang_in(i,j,1,m)*deg2rad ! enddo ! enddo ! enddo call dihang_to_c(dihang_in(1,1,1,1)) ! Store c to cref (they are in COMMON.CHAIN). do k=1,2*nres do kk=1,3 crefjlee(kk,k)=c(kk,k) enddo enddo call contact(.true.,ncont_ref,icont_ref,co) ! do k=1,nres ! write(icsa_native_int,200) (crefjlee(i,k),i=1,3) ! enddo close(icsa_native_int) 200 format (8f10.4) return end subroutine from_int !----------------------------------------------------------------------------- subroutine dihang_to_c(aarray) use geometry_data use csa_data use geometry, only:chainbuild ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.VAR' integer :: i real(kind=8),dimension(mxang,nres,mxch) :: aarray !(mxang,maxres,mxch) ! do i=4,nres ! phi(i)=dihang_in(1,i-2,1,1) ! enddo do i=2,nres-1 theta(i+1)=aarray(1,i,1) phi(i+2)=aarray(2,i,1) alph(i)=aarray(3,i,1) omeg(i)=aarray(4,i,1) enddo call chainbuild return end subroutine dihang_to_c !----------------------------------------------------------------------------- ! geomout.F !----------------------------------------------------------------------------- #ifdef NOXDR subroutine cartout(time) #else subroutine cartoutx(time) #endif use geometry_data, only: c,nres use energy_data use MD_data, only: potE,t_bath ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.IOUNITS' ! include 'COMMON.HEADER' ! include 'COMMON.SBRIDGE' ! include 'COMMON.DISTFIT' ! include 'COMMON.MD' real(kind=8) :: time !el local variables integer :: j,k,i #if defined(AIX) || defined(PGI) open(icart,file=cartname,position="append") #else open(icart,file=cartname,access="append") #endif write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath if (dyn_ss) then write (icart,'(i4,$)') & nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss) else write (icart,'(i4,$)') & nss,(ihpb(j),jhpb(j),j=1,nss) endif write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,& (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),& (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back) write (icart,'(8f10.5)') & ((c(k,j),k=1,3),j=1,nres),& ((c(k,j+nres),k=1,3),j=nnt,nct) close(icart) return #ifdef NOXDR end subroutine cartout #else end subroutine cartoutx #endif !----------------------------------------------------------------------------- #ifndef NOXDR subroutine cartout(time) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use geometry_data, only: c,nres use energy_data use MD_data, only: potE,t_bath #ifdef MPI use MPI_data include 'mpif.h' ! include 'COMMON.SETUP' #else integer,parameter :: me=0 #endif ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.IOUNITS' ! include 'COMMON.HEADER' ! include 'COMMON.SBRIDGE' ! include 'COMMON.DISTFIT' ! include 'COMMON.MD' real(kind=8) :: time integer :: iret,itmp real(kind=4) :: prec real(kind=4),dimension(3,2*nres+2) :: xcoord !(3,maxres2+2) (maxres2=2*maxres !el local variables integer :: j,i,ixdrf #ifdef AIX call xdrfopen_(ixdrf,cartname, "a", iret) call xdrffloat_(ixdrf, real(time), iret) call xdrffloat_(ixdrf, real(potE), iret) call xdrffloat_(ixdrf, real(uconst), iret) call xdrffloat_(ixdrf, real(uconst_back), iret) call xdrffloat_(ixdrf, real(t_bath), iret) call xdrfint_(ixdrf, nss, iret) do j=1,nss if (dyn_ss) then call xdrfint_(ixdrf, idssb(j)+nres, iret) call xdrfint_(ixdrf, jdssb(j)+nres, iret) else call xdrfint_(ixdrf, ihpb(j), iret) call xdrfint_(ixdrf, jhpb(j), iret) endif enddo call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret) do i=1,nfrag call xdrffloat_(ixdrf, real(qfrag(i)), iret) enddo do i=1,npair call xdrffloat_(ixdrf, real(qpair(i)), iret) enddo do i=1,nfrag_back call xdrffloat_(ixdrf, real(utheta(i)), iret) call xdrffloat_(ixdrf, real(ugamma(i)), iret) call xdrffloat_(ixdrf, real(uscdiff(i)), iret) enddo #else call xdrfopen(ixdrf,cartname, "a", iret) call xdrffloat(ixdrf, real(time), iret) call xdrffloat(ixdrf, real(potE), iret) call xdrffloat(ixdrf, real(uconst), iret) call xdrffloat(ixdrf, real(uconst_back), iret) call xdrffloat(ixdrf, real(t_bath), iret) call xdrfint(ixdrf, nss, iret) do j=1,nss if (dyn_ss) then call xdrfint(ixdrf, idssb(j)+nres, iret) call xdrfint(ixdrf, jdssb(j)+nres, iret) else call xdrfint(ixdrf, ihpb(j), iret) call xdrfint(ixdrf, jhpb(j), iret) endif enddo call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret) do i=1,nfrag call xdrffloat(ixdrf, real(qfrag(i)), iret) enddo do i=1,npair call xdrffloat(ixdrf, real(qpair(i)), iret) enddo do i=1,nfrag_back call xdrffloat(ixdrf, real(utheta(i)), iret) call xdrffloat(ixdrf, real(ugamma(i)), iret) call xdrffloat(ixdrf, real(uscdiff(i)), iret) enddo #endif prec=10000.0 do i=1,nres do j=1,3 xcoord(j,i)=c(j,i) enddo enddo do i=nnt,nct do j=1,3 xcoord(j,nres+i-nnt+1)=c(j,i+nres) enddo enddo itmp=nres+nct-nnt+1 #ifdef AIX call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret) call xdrfclose_(ixdrf, iret) #else call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret) call xdrfclose(ixdrf, iret) #endif return end subroutine cartout #endif !----------------------------------------------------------------------------- subroutine statout(itime) use energy_data use control_data use MD_data use MPI_data use compare, only:rms_nac_nnc ! 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.DISTFIT' ! include 'COMMON.MD' ! include 'COMMON.REMD' ! include 'COMMON.SETUP' integer :: itime real(kind=8),dimension(0:n_ene) :: energia ! double precision gyrate !el external gyrate !el common /gucio/ cm character(len=256) :: line1,line2 character(len=4) :: format1,format2 character(len=30) :: format !el local variables integer :: i,ii1,ii2,j real(kind=8) :: rms,frac,frac_nn,co,distance #ifdef AIX if(itime.eq.0) then open(istat,file=statname,position="append") endif #else #ifdef PGI open(istat,file=statname,position="append") #else open(istat,file=statname,access="append") #endif #endif if (AFMlog.gt.0) then if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.false.) write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')& itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),& rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& potEcomp(23),me format1="a133" else !C print *,'A CHUJ',potEcomp(23) write (line1,'(i10,f15.2,7f12.3,i5,$)') & itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),& kinetic_T,t_bath,gyrate(),& potEcomp(23),me format1="a114" endif else if (selfguide.gt.0) then distance=0.0 do j=1,3 distance=distance+(c(j,afmend)-c(j,afmbeg))**2 enddo distance=dsqrt(distance) if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.false.) write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, & f9.3,i5,$)') & itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),& rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& distance,potEcomp(23),me format1="a133" !C print *,"CHUJOWO" else !C print *,'A CHUJ',potEcomp(23) write (line1,'(i10,f15.2,8f12.3,i5,$)')& itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), & kinetic_T,t_bath,gyrate(),& distance,potEcomp(23),me format1="a114" endif else if (velnanoconst.ne.0 ) then if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.false.) write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, & f9.3,i5,$)') & itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),& rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& vecsim,vectrue,me format1="a133" !C print *,"CHUJOWO" else !C print *,'A CHUJ',potEcomp(23) write (line1,'(i10,f15.2,8f12.3,i5,$)')& itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), & kinetic_T,t_bath,gyrate(),& vecsim,vectrue,me format1="a114" endif else if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.false.) write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') & itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),& rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me format1="a133" else write (line1,'(i10,f15.2,7f12.3,i5,$)') & itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),& amax,kinetic_T,t_bath,gyrate(),me format1="a114" endif ENDIF if(usampl.and.totT.gt.eq_time) then write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,& (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),& (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back) write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair & +21*nfrag_back else format2="a001" line2=' ' endif if (print_compon) then if(itime.eq.0) then write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,& ",20a12)" write (istat,format) "#","",& (ename(print_order(i)),i=1,nprint_ene) endif write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,& ",20f12.3)" write (istat,format) line1,line2,& (potEcomp(print_order(i)),i=1,nprint_ene) else write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")" write (istat,format) line1,line2 endif #if defined(AIX) call flush(istat) #else close(istat) #endif return end subroutine statout !----------------------------------------------------------------------------- ! readrtns_CSA.F !----------------------------------------------------------------------------- subroutine readrtns use control_data use energy_data use MPI_data use muca_md, only:read_muca ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.SETUP' ! include 'COMMON.CONTROL' ! include 'COMMON.SBRIDGE' ! include 'COMMON.IOUNITS' logical :: file_exist integer :: i ! Read force-field parameters except weights ! call parmread ! Read job setup parameters call read_control ! Read force-field parameters except weights call parmread ! Read control parameters for energy minimzation if required if (minim) call read_minim ! Read MCM control parameters if required if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread ! Read MD control parameters if reqjuired if (modecalc.eq.12) call read_MDpar ! Read MREMD control parameters if required if (modecalc.eq.14) then call read_MDpar call read_REMDpar endif ! Read MUCA control parameters if required if (lmuca) call read_muca ! Read CSA control parameters if required (from fort.40 if exists ! otherwise from general input file) if (modecalc.eq.8) then inquire (file="fort.40",exist=file_exist) if (.not.file_exist) call csaread endif !fmc if (modecalc.eq.10) call mcmfread ! Read molecule information, molecule geometry, energy-term weights, and ! restraints if requested call molread ! Print restraint information #ifdef MPI if (.not. out1file .or. me.eq.king) then #endif if (nhpb.gt.nss) & write (iout,'(a,i5,a)') "The following",nhpb-nss,& " distance constraints have been imposed" do i=nss+1,nhpb write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i) enddo #ifdef MPI endif #endif ! print *,"Processor",myrank," leaves READRTNS" ! write(iout,*) "end readrtns" return end subroutine readrtns !----------------------------------------------------------------------------- subroutine molread ! ! Read molecular data. ! ! use control, only: ilen use control_data use geometry_data use energy_data use energy use compare_data use MD_data, only: t_bath use MPI_data use compare, only:seq_comp,contact use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer :: error_msg,ierror,ierr,ierrcode #endif ! 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.DBASE' ! include 'COMMON.THREAD' ! include 'COMMON.CONTACTS' ! include 'COMMON.TORCNSTR' ! include 'COMMON.TIME1' ! include 'COMMON.BOUNDS' ! include 'COMMON.MD' ! include 'COMMON.SETUP' character(len=4),dimension(:,:),allocatable :: sequence !(maxres,maxmolecules) ! integer :: rescode ! double precision x(maxvar) character(len=256) :: pdbfile character(len=800) :: weightcard character(len=80) :: weightcard_t!,ucase ! integer,dimension(:),allocatable :: itype_pdb !(maxres) ! common /pizda/ itype_pdb logical :: fail !seq_comp, real(kind=8) :: energia(0:n_ene) ! integer ilen !el external ilen !el local varables integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2,mnum real(kind=8),dimension(3,maxres2+2) :: c_alloc real(kind=8),dimension(3,0:maxres2) :: dc_alloc real(kind=8),dimension(:,:), allocatable :: secprob integer,dimension(maxres) :: itype_alloc integer :: iti,nsi,maxsi,itrial,itmp real(kind=8) :: wlong,scalscp,co,ssscale,phihel,phibet,sigmahel,& sigmabet,sumv,nres_temp allocate(weights(n_ene)) !----------------------------- allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres allocate(dc(3,0:2*maxres)) !(3,0:maxres2) allocate(itype(maxres,5)) !(maxres) allocate(istype(maxres)) ! ! Zero out tables. ! c(:,:)=0.0D0 dc(:,:)=0.0D0 itype(:,:)=0 !----------------------------- ! ! Body ! ! Read weights of the subsequent energy terms. call card_concat(weightcard,.true.) call reada(weightcard,'WLONG',wlong,1.0D0) call reada(weightcard,'WSC',wsc,wlong) call reada(weightcard,'WSCP',wscp,wlong) 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,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) 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) 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,'WSHIELD',wshield,0.05D0) call reada(weightcard,'LIPSCALE',lipscale,1.0D0) call reada(weightcard,'WLT',wliptran,1.0D0) call reada(weightcard,'WTUBE',wtube,1.0d0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) call reada(weightcard,'WSCBASE',wscbase,0.0D0) if (index(weightcard,'SOFT').gt.0) ipot=6 call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0) call reada(weightcard,'WCATCAT',wcatcat,0.0d0) call reada(weightcard,'WCATPROT',wcatprot,0.0d0) call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0) call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0) call reada(weightcard,'WPEPBASE',wpepbase,1.0d0) call reada(weightcard,'WSCPHO',wscpho,0.0d0) call reada(weightcard,'WPEPPHO',wpeppho,0.0d0) ! 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 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)=wpepbase weights(48)=wscpho weights(49)=wpeppho weights(50)=wcatnucl weights(56)=wcat_tran if(me.eq.king.or..not.out1file) & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,& wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,& wturn4,wturn6 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)'/ & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ & 'WTURN6= ',f10.6,' (turns, 6th order)') if(me.eq.king.or..not.out1file)then 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 endif r0_corr=cutoff_corr-delt_corr do i=1,ntyp 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 rescale_weights(t_bath) if(me.eq.king.or..not.out1file) & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,& wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,& wturn4,wturn6 22 format (/'Energy-term weights (scaled):'// & '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)'/ & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ & 'WTURN6= ',f10.6,' (turns, 6th order)') if(me.eq.king.or..not.out1file) & write (iout,*) "Reference temperature for weights calculation:",& temp0 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(me.eq.king.or..not.out1file) then write (iout,*) "Parameters of the SS-bond potential:" write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,& " AKCT",akct write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth write (iout,*)" HT",Ht print *,'indpdb=',indpdb,' pdbref=',pdbref endif if (indpdb.gt.0 .or. pdbref) then read(inp,'(a)') pdbfile if(me.eq.king.or..not.out1file) & write (iout,'(2a)') 'PDB data will be read from file ',& pdbfile(:ilen(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 if (.not.allocated(crefjlee)) allocate (crefjlee(3,2*nres+2)) do i=1,2*nres do j=1,3 crefjlee(j,i)=c(j,i) enddo enddo #ifdef DEBUG do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), & (crefjlee(j,i+nres),j=1,3) enddo #endif ! call int_from_cart1(.true.) ! print *,'Finished reading pdb data' if(me.eq.king.or..not.out1file) & write (iout,'(a,i3,a,i3)')'nsup=',nsup,& ' nstart_sup=',nstart_sup !,"ergwergewrgae" !el if(.not.allocated(itype_pdb)) allocate(itype_pdb(nres)) do i=1,nres itype_pdb(i)=itype(i,1) enddo close (ipdbin) nnt=nstart_sup nct=nstart_sup+nsup-1 !el if(.not.allocated(icont_ref)) allocate(icont_ref(2,(nres/2)*nres)) ! maxcont=12*maxres call contact(.false.,ncont_ref,icont_ref,co) if (sideadd) then if(me.eq.king.or..not.out1file) & write(iout,*)'Adding sidechains' maxsi=1000 do i=2,nres-1 iti=itype(i,1) if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i)) nsi=nsi+1 enddo if(fail) write(iout,*)'Adding sidechain failed for res ',& i,' after ',nsi,' trials' endif enddo endif endif print *,"CZY TU DOCHODZE" if (indpdb.eq.0) then nres_molec(:)=0 allocate(sequence(maxres,5)) ! itype(:,:)=0 itmp=0 if (protein) then ! Read sequence if not taken from the pdb file. molec=1 read (inp,*) nres_molec(molec) print *,'nres=',nres if (iscode.gt.0) then read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec)) else read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec)) endif ! read(inp,*) weightcard_t ! print *,"po seq" weightcard_t ! Convert sequence to numeric code do i=1,nres_molec(molec) itmp=itmp+1 itype(i,1)=rescode(i,sequence(i,molec),iscode,molec) print *,itype(i,1) enddo endif ! read(inp,*) weightcard_t ! print *,"po seq", weightcard_t if (nucleic) then ! Read sequence if not taken from the pdb file. molec=2 read (inp,*) nres_molec(molec) print *,'nres=',nres_molec(molec) ! allocate(sequence(maxres,5)) ! if (iscode.gt.0) then read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec)) print *,"KUR**" print *,(sequence(i,molec),i=1,nres_molec(molec)) ! Convert sequence to numeric code do i=1,nres_molec(molec) itmp=itmp+1 istype(itmp)=sugarcode(sequence(i,molec)(1:1),i) sequence(i,molec)=sequence(i,molec)(1:2) itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec) write(iout,*) "NUCLE=", itype(itmp,molec) enddo endif if (ions) then ! Read sequence if not taken from the pdb file. molec=5 read (inp,*) nres_molec(molec) ! print *,'nres=',nres read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec)) ! Convert sequence to numeric code print *,nres_molec(molec) do i=1,nres_molec(molec) itmp=itmp+1 print *,itmp,"itmp" itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec) enddo endif nres=0 do i=1,5 nres=nres+nres_molec(i) print *,"nres_molec",nres,nres_molec(i) enddo ! Assign initial virtual bond lengths if(.not.allocated(molnum)) then allocate(molnum(nres+1)) itmp=0 do i=1,5 do j=1,nres_molec(i) itmp=itmp+1 molnum(itmp)=i enddo enddo ! print *,nres_molec(i) endif print *,nres,"nres" if(.not.allocated(vbld)) then print *, "I DO ENTER" allocate(vbld(2*nres)) endif if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres)) do i=2,nres mnum=molnum(i) if (molnum(i).eq.1) then vbld(i)=vbl vbld_inv(i)=vblinv else write(iout,*) "typy",itype(i,mnum),ntyp1_molec(mnum),i vbld(i)=6.0 vbld_inv(i)=1.0/6.0 if ((itype(i,mnum).eq.ntyp1_molec(mnum)).or.& (itype(i-1,mnum).eq.ntyp1_molec(mnum))) then vbld(i)=1.9 vbld_inv(i)=1.0/1.9 endif endif enddo do i=2,nres-1 mnum=molnum(i) if (molnum(i).eq.1) then ! print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i vbld(i+nres)=dsc(iabs(itype(i,molnum(i)))) vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i)))) else write(iout,*) "typy2",itype(i,mnum),ntyp1_molec(mnum),i if (itype(i,mnum).eq.ntyp1_molec(mnum)) then vbld_inv(i+nres)=1.0 vbld(i+nres)=0.0 else vbld(i+nres)=vbldsc0_nucl(1,iabs(itype(i,molnum(i)))) vbld_inv(i+nres)=1.0/vbldsc0_nucl(1,iabs(itype(i,molnum(i)))) endif endif ! write (iout,*) "i",i," itype",itype(i,1), ! & " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres) enddo endif ! print *,nres ! print '(20i4)',(itype(i,1),i=1,nres) !---------------------------- !el reallocate tables ! do i=1,maxres2 ! do j=1,3 ! c_alloc(j,i)=c(j,i) ! dc_alloc(j,i)=dc(j,i) ! enddo ! enddo ! do i=1,maxres !elwrite(iout,*) "itype",i,itype(i,1) ! itype_alloc(i)=itype(i,1) ! enddo ! deallocate(c) ! deallocate(dc) ! deallocate(itype) ! allocate(c(3,2*nres+4)) ! allocate(dc(3,0:2*nres+2)) ! allocate(itype(nres+2)) allocate(itel(nres+2)) itel(:)=0 ! do i=1,2*nres+2 ! do j=1,3 ! c(j,i)=c_alloc(j,i) ! dc(j,i)=dc_alloc(j,i) ! enddo ! enddo ! do i=1,nres+2 ! itype(i,1)=itype_alloc(i) ! itel(i)=0 ! enddo !-------------------------- do i=1,nres #ifdef PROCOR if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then #else if (itype(i,1).eq.ntyp1) 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 if(me.eq.king.or..not.out1file)then write (iout,*) "ITEL" print *,nres,"nres" do i=1,nres-1 write (iout,*) i,itype(i,1),itel(i) enddo print *,'Call Read_Bridge.' endif call read_bridge !-------------------------------- ! print *,"tu dochodze" ! znamy nres oraz nss można zaalokowac potrzebne tablice call alloc_geo_arrays call alloc_ener_arrays !-------------------------------- ! 8/13/98 Set limits to generating the dihedral angles do i=1,nres phibound(1,i)=-pi phibound(2,i)=pi enddo read (inp,*) ndih_constr if (ndih_constr.gt.0) then raw_psipred=.false. allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr) allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr) allocate(ftors(ndih_constr)) !(maxdih_constr) ! read (inp,*) ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), & i=1,ndih_constr) if(me.eq.king.or..not.out1file)then write (iout,*) & 'There are',ndih_constr,' constraints on phi angles.' do i=1,ndih_constr write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), & ftors(i) enddo endif do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo ! if(me.eq.king.or..not.out1file) & ! write (iout,*) 'FTORS',ftors do i=1,ndih_constr ii = idih_constr(i) phibound(1,ii) = phi0(i)-drange(i) phibound(2,ii) = phi0(i)+drange(i) enddo else if (ndih_constr.lt.0) then raw_psipred=.true. allocate(idih_constr(nres)) allocate(secprob(3,nres)) allocate(vpsipred(3,nres)) allocate(sdihed(2,nres)) call card_concat(weightcard,.true.) call reada(weightcard,"PHIHEL",phihel,50.0D0) call reada(weightcard,"PHIBET",phibet,180.0D0) call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0) call reada(weightcard,"SIGMABET",sigmabet,40.0d0) call reada(weightcard,"WDIHC",wdihc,0.591D0) write (iout,*) "Weight of dihedral angle restraints",wdihc read(inp,'(9x,3f7.3)') & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct) write (iout,*) "The secprob array" do i=nnt,nct write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3) enddo ndih_constr=0 do i=nnt+3,nct if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 & .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then ndih_constr=ndih_constr+1 idih_constr(ndih_constr)=i sumv=0.0d0 do j=1,3 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2) sumv=sumv+vpsipred(j,ndih_constr) enddo do j=1,3 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv enddo phibound(1,ndih_constr)=phihel*deg2rad phibound(2,ndih_constr)=phibet*deg2rad sdihed(1,ndih_constr)=sigmahel*deg2rad sdihed(2,ndih_constr)=sigmabet*deg2rad endif enddo endif if (with_theta_constr) then !C with_theta_constr is keyword allowing for occurance of theta constrains read (inp,*) ntheta_constr !C ntheta_constr is the number of theta constrains if (ntheta_constr.gt.0) then !C read (inp,*) ftors allocate(itheta_constr(ntheta_constr)) allocate(theta_constr0(ntheta_constr)) allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr)) read (inp,*) (itheta_constr(i),theta_constr0(i), & theta_drange(i),for_thet_constr(i), & i=1,ntheta_constr) !C the above code reads from 1 to ntheta_constr !C itheta_constr(i) residue i for which is theta_constr !C theta_constr0 the global minimum value !C theta_drange is range for which there is no energy penalty !C for_thet_constr is the force constant for quartic energy penalty !C E=k*x**4 if(me.eq.king.or..not.out1file)then write (iout,*) & 'There are',ntheta_constr,' constraints on phi angles.' do i=1,ntheta_constr write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), & theta_drange(i), & for_thet_constr(i) enddo endif do i=1,ntheta_constr theta_constr0(i)=deg2rad*theta_constr0(i) theta_drange(i)=deg2rad*theta_drange(i) enddo !C if(me.eq.king.or..not.out1file) !C & write (iout,*) 'FTORS',ftors !C do i=1,ntheta_constr !C ii = itheta_constr(i) !C thetabound(1,ii) = phi0(i)-drange(i) !C thetabound(2,ii) = phi0(i)+drange(i) !C enddo endif ! ntheta_constr.gt.0 endif! with_theta_constr nnt=1 #ifdef MPI if (me.eq.king) then #endif write (iout,'(a)') 'Boundaries in phi angle sampling:' do i=1,nres write (iout,'(a3,i5,2f10.1)') & restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg enddo #ifdef MP endif #endif 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 (pdbref) then if(me.eq.king.or..not.out1file) & 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,1),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,1),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 if (nsup.eq.0) nsup=nct-nnt+1 if (nstart_sup.eq.0) nstart_sup=nnt if (nstart_seq.eq.0) nstart_seq=nnt if(me.eq.king.or..not.out1file) & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,& ' nstart_seq=',nstart_seq !,"242343453254" endif !--- Zscore rms ------- if (nz_start.eq.0) nz_start=nnt if (nz_end.eq.0 .and. nsup.gt.0) then nz_end=nnt+nsup-1 else if (nz_end.eq.0) then nz_end=nct endif if(me.eq.king.or..not.out1file)then write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end write (iout,*) 'IZ_SC=',iz_sc, 'NCT=',nct endif !---------------------- call init_int_table if (refstr) then if (.not.pdbref) then call read_angles(inp,*38) goto 39 38 write (iout,'(a)') 'Error reading reference structure.' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,IERROR) stop 'Error reading reference structure' #endif 39 call chainbuild call setup_var !zscore call geom_to_var(nvar,coord_exp_zs(1,1)) 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 call contact(.true.,ncont_ref,icont_ref,co) endif ! write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup ! call flush(iout) !EL if (constr_dist.gt.0) call read_dist_constr !EL write (iout,*) "After read_dist_constr nhpb",nhpb !EL if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp !EL call hpb_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co if (pdbref) then if(me.eq.king.or..not.out1file) & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup do i=1,ncont_ref do j=1,2 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup enddo if(me.eq.king.or..not.out1file) & write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',& icont_ref(1,i),' ',& restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i) enddo endif if (constr_homology.gt.0) then ! write (iout,*) "Calling read_constr_homology" ! call flush(iout) call read_constr_homology if (indpdb.gt.0 .or. pdbref) then do i=1,2*nres do j=1,3 c(j,i)=crefjlee(j,i) cref(j,i,1)=crefjlee(j,i) enddo enddo endif #define DEBUG #ifdef DEBUG write (iout,*) "sc_loc_geom: Array C" do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),& (c(j,i+nres),j=1,3) enddo write (iout,*) "Array Cref" do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),& (cref(j,i+nres,1),j=1,3) enddo #endif call int_from_cart1(.false.) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) enddo enddo do i=2,nres-1 do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) enddo enddo else homol_nset=0 if (start_from_model) then nmodel_start=0 do read(inp,'(a)',end=332,err=332) pdbfile if (me.eq.king .or. .not. out1file)& write (iout,'(a,5x,a)') 'Opening PDB file',& pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=336) goto 335 336 write (iout,'(a,5x,a)') 'Error opening PDB file',& pdbfile(:ilen(pdbfile)) call flush(iout) stop 335 continue unres_pdb=.false. nres_temp=nres ! call readpdb call readpdb_template(nmodel_start+1) close(ipdbin) if (nres.ge.nres_temp) then nmodel_start=nmodel_start+1 pdbfiles_chomo(nmodel_start)=pdbfile do i=1,2*nres do j=1,3 chomo(j,i,nmodel_start)=c(j,i) enddo enddo else if (me.eq.king .or. .not. out1file) & write (iout,'(a,2i7,1x,a)') & "Different number of residues",nres_temp,nres, & " model skipped." endif nres=nres_temp enddo 332 continue if (nmodel_start.eq.0) then if (me.eq.king .or. .not. out1file) & write (iout,'(a)') & "No valid starting model found START_FROM_MODELS is OFF" start_from_model=.false. endif write (iout,*) "nmodel_start",nmodel_start endif endif endif if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp if (velnanoconst.ne.0) call read_afmnano call hpb_partition if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. & modecalc.ne.10) then ! If input structure hasn't been supplied from the PDB file read or generate ! initial geometry. if (iranconf.eq.0 .and. .not. extconf) then if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) & write (iout,'(a)') 'Initial geometry will be read in.' if (read_cart) then read(inp,'(8f10.5)',end=36,err=36) & ((c(l,k),l=1,3),k=1,nres),& ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "Exit READ_CART" 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) call int_from_cart1(.true.) write (iout,*) "Finish INT_TO_CART" do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1) enddo enddo do i=nnt,nct if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres) enddo endif enddo return else write(iout,*) "read angles from input" call read_angles(inp,*36) call chainbuild endif goto 37 36 write (iout,'(a)') 'Error reading angle file.' #ifdef MPI call mpi_finalize( MPI_COMM_WORLD,IERR ) #endif stop 'Error reading angle file.' 37 continue else if (extconf) then if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) & write (iout,'(a)') 'Extended chain initial geometry.' do i=3,nres theta(i)=90d0*deg2rad if (molnum(i).eq.2) theta(i)=160d0*deg2rad enddo do i=4,nres phi(i)=180d0*deg2rad if (molnum(i).eq.2) phi(i)=30d0*deg2rad enddo do i=2,nres-1 alph(i)=110d0*deg2rad if (molnum(i).eq.2) alph(i)=30d0*deg2rad enddo do i=2,nres-1 omeg(i)=-120d0*deg2rad if (molnum(i).eq.2) omeg(i)=60d0*deg2rad if (itype(i,1).le.0) omeg(i)=-omeg(i) enddo call chainbuild else if(me.eq.king.or..not.out1file) & write (iout,'(a)') 'Random-generated initial geometry.' #ifdef MPI if (me.eq.king .or. fg_rank.eq.0 .and. & ( modecalc.eq.12 .or. modecalc.eq.14) ) then #endif do itrial=1,100 itmp=1 call gen_rand_conf(itmp,*30) goto 40 30 write (iout,*) 'Failed to generate random conformation',& ', itrial=',itrial write (*,*) 'Processor:',me,& ' Failed to generate random conformation',& ' itrial=',itrial call intout #ifdef AIX call flush_(iout) #else call flush(iout) #endif enddo write (iout,'(a,i3,a)') 'Processor:',me,& ' error in generating random conformation.' write (*,'(a,i3,a)') 'Processor:',me,& ' error in generating random conformation.' call flush(iout) #ifdef MPI call MPI_Abort(mpi_comm_world,error_msg,ierrcode) 40 continue endif #else do itrial=1,100 itmp=1 call gen_rand_conf(itmp,*335) goto 40 335 write (iout,*) 'Failed to generate random conformation',& ', itrial=',itrial write (*,*) 'Failed to generate random conformation',& ', itrial=',itrial enddo write (iout,'(a,i3,a)') 'Processor:',me,& ' error in generating random conformation.' write (*,'(a,i3,a)') 'Processor:',me,& ' error in generating random conformation.' stop 40 continue #endif endif elseif (modecalc.eq.4) then read (inp,'(a)') intinname open (intin,file=intinname,status='old',err=333) if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) & write (iout,'(a)') 'intinname',intinname write (*,'(a)') 'Processor',myrank,' intinname',intinname goto 334 333 write (iout,'(2a)') 'Error opening angle file ',intinname #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,IERR) #endif stop 'Error opening angle file.' 334 continue endif ! Generate distance constraints, if the PDB structure is to be regularized. if (nthread.gt.0) then call read_threadbase endif call setup_var if (me.eq.king .or. .not. out1file) & call intout if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then write (iout,'(/a,i3,a)') & 'The chain contains',ns,' disulfide-bridging cysteines.' write (iout,'(20i4)') (iss(i),i=1,ns) if (dyn_ss) then write(iout,*)"Running with dynamic disulfide-bond formation" else write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss i1=ihpb(i)-nres i2=jhpb(i)-nres it1=itype(i1,1) it2=itype(i2,1) if (me.eq.king.or..not.out1file) & write (iout,'(2a,i3,3a,i3,a,3f10.3)') & restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),& ebr,forcon(i) enddo write (iout,'(a)') endif endif 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 call hpb_partition do i=1,ns dyn_ss_mask(iss(i))=.true. enddo endif if (i2ndstr.gt.0) call secstrp2dihc if (indpdb.gt.0) then write(iout,*) "WCHODZE TU!!" call int_from_cart1(.true.) endif ! call geom_to_var(nvar,x) ! call etotal(energia(0)) ! call enerprint(energia(0)) ! call briefout(0,etot) ! stop !d write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT !d write (iout,'(a)') 'Variable list:' !d write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar) #ifdef MPI if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') & 'Processor',myrank,': end reading molecular data.' #endif return end subroutine molread !----------------------------------------------------------------------------- subroutine read_constr_homology use control, only:init_int_table,homology_partition use MD_data, only:iset ! implicit none ! include 'DIMENSIONS' !#ifdef MPI ! include 'mpif.h' !#endif ! include 'COMMON.SETUP' ! include 'COMMON.CONTROL' ! include 'COMMON.HOMOLOGY' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.MD' ! include 'COMMON.QRESTR' ! include 'COMMON.GEO' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.VAR' ! ! double precision odl_temp,sigma_odl_temp,waga_theta,waga_d, ! & dist_cut ! common /przechowalnia/ odl_temp(maxres,maxres,max_template), ! & sigma_odl_temp(maxres,maxres,max_template) character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard integer :: ki,i,ii,j,k,l integer, dimension (:), allocatable :: ii_in_use integer :: i_tmp,idomain_tmp,& irec,ik,iistart,nres_temp ! integer :: iset ! external :: ilen logical :: liiflag,lfirst integer :: i01,i10 ! ! FP - Nov. 2014 Temporary specifications for new vars ! real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,& rescore3_tmp, dist_cut real(kind=8), dimension (:,:),allocatable :: rescore real(kind=8), dimension (:,:),allocatable :: rescore2 real(kind=8), dimension (:,:),allocatable :: rescore3 real(kind=8) :: distal character*24 tpl_k_rescore character*256 pdbfile ! ----------------------------------------------------------------- ! Reading multiple PDB ref structures and calculation of retraints ! not using pre-computed ones stored in files model_ki_{dist,angle} ! FP (Nov., 2014) ! ----------------------------------------------------------------- ! ! ! Alternative: reading from input call card_concat(controlcard,.true.) call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0) call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0) call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) call readi(controlcard,"HOMOL_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0) if(.not.read2sigma.and.start_from_model) then if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)& write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA' start_from_model=.false. iranconf=(indpdb.le.0) endif if(start_from_model .and. (me.eq.king .or. .not. out1file))& write(iout,*) 'START_FROM_MODELS is ON' ! if(start_from_model .and. rest) then ! if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then ! write(iout,*) 'START_FROM_MODELS is OFF' ! write(iout,*) 'remove restart keyword from input' ! endif ! endif if (start_from_model) nmodel_start=constr_homology if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset)) if (homol_nset.gt.1)then call card_concat(controlcard,.true.) read(controlcard,*) (waga_homology(i),i=1,homol_nset) if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then ! write(iout,*) "iset homology_weight " do i=1,homol_nset write(iout,*) i,waga_homology(i) enddo endif iset=mod(kolor,homol_nset)+1 else iset=1 waga_homology(1)=1.0 endif !d write (iout,*) "nnt",nnt," nct",nct !d call flush(iout) if (read_homol_frag) then call read_klapaucjusz else lim_odl=0 lim_dih=0 ! ! write(iout,*) 'nnt=',nnt,'nct=',nct ! ! do i = nnt,nct ! do k=1,constr_homology ! idomain(k,i)=0 ! enddo ! enddo idomain=0 ! ii=0 ! do i = nnt,nct-2 ! do j=i+2,nct ! ii=ii+1 ! ii_in_use(ii)=0 ! enddo ! enddo ii_in_use=0 if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology)) if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology)) do k=1,constr_homology read(inp,'(a)') pdbfile pdbfiles_chomo(k)=pdbfile if(me.eq.king .or. .not. out1file) & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',& pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file',& pdbfile(:ilen(pdbfile)) stop 34 continue ! print *,'Begin reading pdb data' ! ! Files containing res sim or local scores (former containing sigmas) ! write(kic2,'(bz,i2.2)') k tpl_k_rescore="template"//kic2//".sco" write(iout,*) "tpl_k_rescore=",tpl_k_rescore unres_pdb=.false. nres_temp=nres write(iout,*) "read2sigma",read2sigma if (read2sigma) then call readpdb_template(k) else call readpdb endif write(iout,*) "after readpdb" if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology)) nres_chomo(k)=nres nres=nres_temp if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres)) if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres)) if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres)) if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1))) if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres)) if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres)) if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200)) if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200)) if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200)) if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200)) if(.not.allocated(dih)) allocate(dih(constr_homology,nres)) if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres)) if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres)) if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres)) ! if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres)) if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres)) if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres)) if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres)) if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres)) ! if(.not.allocated(distance)) allocate(distance(constr_homology)) ! if(.not.allocated(distancek)) allocate(distancek(constr_homology)) ! ! Distance restraints ! ! ... --> odl(k,ii) ! Copy the coordinates from reference coordinates (?) do i=1,2*nres_chomo(k) do j=1,3 c(j,i)=cref(j,i,1) ! write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo ! ! From read_dist_constr (commented out 25/11/2014 <-> res sim) ! ! write(iout,*) "tpl_k_rescore - ",tpl_k_rescore open (ientin,file=tpl_k_rescore,status='old') if (nnt.gt.1) rescore(k,1)=0.0d0 do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,& rescore3_tmp,idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp rescore3(k,i_tmp)=rescore3_tmp if (.not. out1file .or. me.eq.king)& write(iout,'(a7,i5,3f10.5,i5)') "rescore",& i_tmp,rescore2_tmp,rescore_tmp,& rescore3_tmp,idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp ! rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores ! write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) endif enddo 1401 continue close (ientin) if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 do j=i+2,nct x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) ! write (iout,*) k,i,j,distal,dist2_cut if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. ! write (iout,*) "k",k ! write (iout,*) "i",i," j",j," constr_homology", ! & constr_homology ires_homo(ii)=i jres_homo(ii)=j odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 do ik=i,j sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) enddo sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ! ii=ii+1 ! l_homo(k,ii)=.false. endif enddo enddo lim_odl=ii endif ! write (iout,*) "Distance restraints set" ! call flush(iout) ! ! Theta, dihedral and SC retraints ! if (waga_angle.gt.0.0d0) then ! open (ientin,file=tpl_k_sigma_dih,status='old') ! do irec=1,maxres-3 ! loop for reading sigma_dih ! read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for? ! if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right? ! sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity ! & sigma_dih(k,i+nnt-1) ! enddo !1402 continue ! close (ientin) do i = nnt+3,nct if (idomain(k,i).eq.0) then sigma_dih(k,i)=0.0 cycle endif dih(k,i)=phiref(i) ! right? ! read (ientin,*) sigma_dih(k,i) ! original variant ! write (iout,*) "dih(",k,i,") =",dih(k,i) ! write(iout,*) "rescore(",k,i,") =",rescore(k,i), ! & "rescore(",k,i-1,") =",rescore(k,i-1), ! & "rescore(",k,i-2,") =",rescore(k,i-2), ! & "rescore(",k,i-3,") =",rescore(k,i-3) sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 ! if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 ! write (iout,*) "Raw sigmas for dihedral angle restraints" ! write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i) ! sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* ! rescore(k,i-2)*rescore(k,i-3) ! right expression ? ! Instead of res sim other local measure of b/b str reliability possible if (sigma_dih(k,i).ne.0) & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) ! sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) enddo lim_dih=nct-nnt-2 endif ! write (iout,*) "Dihedral angle restraints set" ! call flush(iout) if (waga_theta.gt.0.0d0) then ! open (ientin,file=tpl_k_sigma_theta,status='old') ! do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds? ! read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for? ! sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity ! & sigma_theta(k,i+nnt-1) ! enddo !1403 continue ! close (ientin) do i = nnt+2,nct ! right? without parallel. ! do i = i=1,nres ! alternative for bounds acc to readpdb? ! do i=ithet_start,ithet_end ! with FG parallel. if (idomain(k,i).eq.0) then sigma_theta(k,i)=0.0 cycle endif thetatpl(k,i)=thetaref(i) ! write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) ! write(iout,*) "rescore(",k,i,") =",rescore(k,i), ! & "rescore(",k,i-1,") =",rescore(k,i-1), ! & "rescore(",k,i-2,") =",rescore(k,i-2) ! read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 ! if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 if (sigma_theta(k,i).ne.0) & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) ! sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* ! rescore(k,i-2) ! right expression ? ! sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif ! write (iout,*) "Angle restraints set" ! call flush(iout) if (waga_d.gt.0.0d0) then ! open (ientin,file=tpl_k_sigma_d,status='old') ! do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds? ! read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for? ! sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity ! & sigma_d(k,i+nnt-1) ! enddo !1404 continue do i = nnt,nct ! right? without parallel. ! do i=2,nres-1 ! alternative for bounds acc to readpdb? ! do i=loc_start,loc_end ! with FG parallel. if (itype(i,1).eq.10) cycle if (idomain(k,i).eq.0 ) then sigma_d(k,i)=0.0 cycle endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) ! write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i) ! write (iout,*) "yytpl(",k,i,") =",yytpl(k,i) ! write (iout,*) "zztpl(",k,i,") =",zztpl(k,i) ! write(iout,*) "rescore(",k,i,") =",rescore(k,i) sigma_d(k,i)=rescore3(k,i) ! right expression ? if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) ! sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? ! sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) ! read (ientin,*) sigma_d(k,i) ! 1st variant enddo endif enddo ! write (iout,*) "SC restraints set" ! call flush(iout) ! ! remove distance restraints not used in any model from the list ! shift data in all arrays ! ! write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. lfirst=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 ! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 ! & .and. distal.le.dist2_cut ) then ! write (iout,*) "i",i," j",j," ii",ii ! call flush(iout) if (ii_in_use(ii).eq.0.and.liiflag.or. & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then liiflag=.false. i10=ii if (lfirst) then lfirst=.false. iistart=ii else if(i10.eq.lim_odl) i10=i10+1 do ki=0,i10-i01-1 ires_homo(iistart+ki)=ires_homo(ki+i01) jres_homo(iistart+ki)=jres_homo(ki+i01) ii_in_use(iistart+ki)=ii_in_use(ki+i01) do k=1,constr_homology odl(k,iistart+ki)=odl(k,ki+i01) sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01) l_homo(k,iistart+ki)=l_homo(k,ki+i01) enddo enddo iistart=iistart+i10-i01 endif endif if (ii_in_use(ii).ne.0.and..not.liiflag) then i01=ii liiflag=.true. endif enddo enddo lim_odl=iistart-1 endif ! write (iout,*) "Removing distances completed" ! call flush(iout) endif ! .not. klapaucjusz if (constr_homology.gt.0) call homology_partition write (iout,*) "After homology_partition" call flush(iout) if (constr_homology.gt.0) call init_int_table write (iout,*) "After init_int_table" call flush(iout) ! endif ! .not. klapaucjusz ! endif ! if (constr_homology.gt.0) call homology_partition ! write (iout,*) "After homology_partition" ! call flush(iout) ! if (constr_homology.gt.0) call init_int_table ! write (iout,*) "After init_int_table" ! call flush(iout) ! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end ! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end ! ! Print restraints ! if (.not.out_template_restr) return !d write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,*) "Distance restraints from templates" do ii=1,lim_odl write(iout,'(3i7,100(2f8.2,1x,l1,4x))') & ii,ires_homo(ii),jres_homo(ii),& (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),& ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" do i=nnt+3,nct write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),& (rad2deg*dih(ki,i),& rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) enddo write (iout,*) "Virtual-bond angle restraints from templates" do i=nnt+2,nct write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),& (rad2deg*thetatpl(ki,i),& rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) enddo write (iout,*) "SC restraints from templates" do i=nnt,nct write(iout,'(i7,100(4f8.2,4x))') i,& (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) enddo endif return end subroutine read_constr_homology !----------------------------------------------------------------------------- subroutine read_klapaucjusz use energy_data implicit none ! include 'DIMENSIONS' !#ifdef MPI ! include 'mpif.h' !#endif ! include 'COMMON.SETUP' ! include 'COMMON.CONTROL' ! include 'COMMON.HOMOLOGY' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.MD' ! include 'COMMON.GEO' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' character(len=256) fragfile integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,& ii_in_use integer, dimension(:,:), allocatable :: iresclust,inclust integer :: nclust character(len=2) :: kic2 character(len=24) :: model_ki_dist, model_ki_angle character(len=500) :: controlcard integer :: ki, i, j, jj,k, l, i_tmp,& idomain_tmp,& ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,& i01,i10,nnt_chain,nct_chain real(kind=8) :: distal logical :: lprn = .true. integer :: nres_temp ! integer :: ilen ! external :: ilen logical :: liiflag,lfirst real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut real(kind=8), dimension (:,:), allocatable :: rescore real(kind=8), dimension (:,:), allocatable :: rescore2 character(len=24) :: tpl_k_rescore character(len=256) :: pdbfile ! ! For new homol impl ! ! include 'COMMON.VAR' ! ! write (iout,*) "READ_KLAPAUCJUSZ" ! print *,"READ_KLAPAUCJUSZ" ! call flush(iout) call getenv("FRAGFILE",fragfile) write (iout,*) "Opening", fragfile call flush(iout) open(ientin,file=fragfile,status="old",err=10) ! write (iout,*) " opened" ! call flush(iout) sigma_theta=0.0 sigma_d=0.0 sigma_dih=0.0 l_homo = .false. nres_temp=nres itype_temp(:)=itype(:,1) ii=0 lim_odl=0 ! write (iout,*) "Entering loop" ! call flush(iout) DO IGR = 1,NCHAIN_GROUP ! write (iout,*) "igr",igr call flush(iout) read(ientin,*) constr_homology,nclust if (start_from_model) then nmodel_start=constr_homology else nmodel_start=0 endif ii_old=lim_odl ichain=iequiv(1,igr) nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1 nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1 ! write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain ! Read pdb files if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology)) if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology)) do k=1,constr_homology read(ientin,'(a)') pdbfile write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', & pdbfile(:ilen(pdbfile)) pdbfiles_chomo(k)=pdbfile open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file',& pdbfile(:ilen(pdbfile)) stop 34 continue unres_pdb=.false. ! nres_temp=nres call readpdb_template(k) nres_chomo(k)=nres ! nres=nres_temp do i=1,nres rescore(k,i)=0.2d0 rescore2(k,i)=1.0d0 enddo enddo ! Read clusters do i=1,nclust read(ientin,*) ninclust(i),nresclust(i) read(ientin,*) (inclust(k,i),k=1,ninclust(i)) read(ientin,*) (iresclust(k,i),k=1,nresclust(i)) enddo ! ! Loop over clusters ! do l=1,nclust do ll = 1,ninclust(l) k = inclust(ll,l) ! write (iout,*) "l",l," ll",ll," k",k do i=1,nres idomain(k,i)=0 enddo do i=1,nresclust(l) if (nnt.gt.1) then idomain(k,iresclust(i,l)+1) = 1 else idomain(k,iresclust(i,l)) = 1 endif enddo ! ! Distance restraints ! ! ... --> odl(k,ii) ! Copy the coordinates from reference coordinates (?) ! nres_temp=nres nres=nres_chomo(k) do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,k) ! write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo ! nres=nres_temp if (waga_dist.ne.0.0d0) then ii=ii_old ! do i = nnt,nct-2 do i = nnt_chain,nct_chain-2 ! do j=i+2,nct do j=i+2,nct_chain x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) ! write (iout,*) k,i,j,distal,dist2_cut if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. ! write (iout,*) "k",k ! write (iout,*) "i",i," j",j," constr_homology", ! & constr_homology ires_homo(ii)=i+chain_border1(1,igr)-1 jres_homo(ii)=j+chain_border1(1,igr)-1 odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 do ik=i,j sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) enddo sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ii=ii+1 ! l_homo(k,ii)=.false. endif enddo enddo lim_odl=ii endif ! ! Theta, dihedral and SC retraints ! if (waga_angle.gt.0.0d0) then do i = nnt_chain+3,nct_chain iii=i+chain_border1(1,igr)-1 if (idomain(k,i).eq.0) then ! sigma_dih(k,i)=0.0 cycle endif dih(k,iii)=phiref(i) sigma_dih(k,iii)= & (rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 ! write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), ! & " sigma_dihed",sigma_dih(k,i) if (sigma_dih(k,iii).ne.0) & sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii)) enddo ! lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then do i = nnt_chain+2,nct_chain iii=i+chain_border1(1,igr)-1 if (idomain(k,i).eq.0) then ! sigma_theta(k,i)=0.0 cycle endif thetatpl(k,iii)=thetaref(i) sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 if (sigma_theta(k,iii).ne.0) & sigma_theta(k,iii)=1.0d0/ & (sigma_theta(k,iii)*sigma_theta(k,iii)) enddo endif if (waga_d.gt.0.0d0) then do i = nnt_chain,nct_chain iii=i+chain_border1(1,igr)-1 if (itype(i,1).eq.10) cycle if (idomain(k,i).eq.0 ) then ! sigma_d(k,i)=0.0 cycle endif xxtpl(k,iii)=xxref(i) yytpl(k,iii)=yyref(i) zztpl(k,iii)=zzref(i) sigma_d(k,iii)=rescore(k,i) if (sigma_d(k,iii).ne.0) & sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii)) ! if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 enddo endif enddo ! l enddo ! ll ! ! remove distance restraints not used in any model from the list ! shift data in all arrays ! ! write (iout,*) "ii_old",ii_old if (waga_dist.ne.0.0d0) then #ifdef DEBUG write (iout,*) "Distance restraints from templates" do iii=1,lim_odl write(iout,'(4i5,100(2f8.2,1x,l1,4x))') & iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), & (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), & ki=1,constr_homology) enddo #endif ii=ii_old liiflag=.true. lfirst=.true. do i=nnt_chain,nct_chain-2 do j=i+2,nct_chain ii=ii+1 ! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 ! & .and. distal.le.dist2_cut ) then ! write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii) ! call flush(iout) if (ii_in_use(ii).eq.0.and.liiflag.or. & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then liiflag=.false. i10=ii if (lfirst) then lfirst=.false. iistart=ii else if(i10.eq.lim_odl) i10=i10+1 do ki=0,i10-i01-1 ires_homo(iistart+ki)=ires_homo(ki+i01) jres_homo(iistart+ki)=jres_homo(ki+i01) ii_in_use(iistart+ki)=ii_in_use(ki+i01) do k=1,constr_homology odl(k,iistart+ki)=odl(k,ki+i01) sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01) l_homo(k,iistart+ki)=l_homo(k,ki+i01) enddo enddo iistart=iistart+i10-i01 endif endif if (ii_in_use(ii).ne.0.and..not.liiflag) then i01=ii liiflag=.true. endif enddo enddo lim_odl=iistart-1 endif lll=lim_odl-ii_old do i=2,nequiv(igr) ichain=iequiv(i,igr) do j=nnt_chain,nct_chain jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr)) do k=1,constr_homology dih(k,jj)=dih(k,j) sigma_dih(k,jj)=sigma_dih(k,j) thetatpl(k,jj)=thetatpl(k,j) sigma_theta(k,jj)=sigma_theta(k,j) xxtpl(k,jj)=xxtpl(k,j) yytpl(k,jj)=yytpl(k,j) zztpl(k,jj)=zztpl(k,j) sigma_d(k,jj)=sigma_d(k,j) enddo enddo jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr)) ! write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj do j=ii_old+1,lim_odl ires_homo(j+lll)=ires_homo(j)+jj jres_homo(j+lll)=jres_homo(j)+jj do k=1,constr_homology odl(k,j+lll)=odl(k,j) sigma_odl(k,j+lll)=sigma_odl(k,j) l_homo(k,j+lll)=l_homo(k,j) enddo enddo ii_old=ii_old+lll lim_odl=lim_odl+lll enddo ENDDO ! IGR if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2 nres=nres_temp itype(:,1)=itype_temp(:) return 10 stop "Error in fragment file" end subroutine read_klapaucjusz !----------------------------------------------------------------------------- end module io