include 'COMMON.CONTROL'
include 'COMMON.SBRIDGE'
include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
logical file_exist
C Read force-field parameters except weights
call parmread
sideadd=(index(controlcard,'SIDEADD').gt.0)
energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
outpdb=(index(controlcard,'PDBOUT').gt.0)
+ outx=(index(controlcard,'XOUT').gt.0)
outmol2=(index(controlcard,'MOL2OUT').gt.0)
pdbref=(index(controlcard,'PDBREF').gt.0)
refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
large = index(controlcard,"LARGE").gt.0
print_compon = index(controlcard,"PRINT_COMPON").gt.0
rattle = index(controlcard,"RATTLE").gt.0
+ preminim = index(controlcard,"PREMINIM").gt.0
+ if (preminim) then
+ dccart=(index(controlcard,'CART').gt.0)
+ call read_minim
+ endif
c if performing umbrella sampling, fragments constrained are read from the fragment file
nset=0
if(usampl) then
write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
if (rattle) write (iout,'(a60)')
& "Rattle algorithm used to constrain the virtual bonds"
+ if (preminim .or. iranconf.gt.0) then
+ write (iout,'(a60)')
+ & "Initial structure will be energy-minimized"
+ endif
endif
reset_fricmat=1000
if (lang.gt.0) then
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
c print *,'Finished reading pdb data'
if(me.eq.king.or..not.out1file)
& write (iout,'(a,i3,a,i3)')'nsup=',nsup,
C READ fragment information!!
C both routines should be in dfa.F file!!
+#ifdef DFA
if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
& wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
call init_dfa_vars
call read_dfa_info
print*, 'read_dfa_info finished!'
endif
+#endif
C
if (pdbref) then
if(me.eq.king.or..not.out1file)
enddo
enddo
endif
+#ifdef DEBUG
+ write (iout,*) "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),j=1,3),
+ & (cref(j,i+nres),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
endif
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,*) time,potE,uconst,t_bath,
+ & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
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)
enddo
endif
enddo
- return
+c return
else
call read_angles(inp,*36)
endif
if (nthread.gt.0) then
call read_threadbase
endif
+ write (iout,*) "READRTNS: Calling setup_var"
+ call flush(iout)
call setup_var
if (me.eq.king .or. .not. out1file)
& call intout
& //'.pdb'
mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
& liczba(:ilen(liczba))//'.mol2'
+ cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
+ & liczba(:ilen(liczba))//'.x'
statname=prefix(:lenpre)//'_'//pot(:lenpot)//
& liczba(:ilen(liczba))//'.stat'
if (lentmp.gt.0)
intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
+ cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
if (lentmp.gt.0)
& call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
include 'COMMON.IOUNITS'
include 'COMMON.MD'
include 'COMMON.CONTROL'
+ integer iset1
read(inp,*) nset,nfrag,npair,nfrag_back
if(me.eq.king.or..not.out1file)
& write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
& " nfrag_back",nfrag_back
- do iset=1,nset
- read(inp,*) mset(iset)
+ do iset1=1,nset
+ read(inp,*) mset(iset1)
do i=1,nfrag
- read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
- & qinfrag(i,iset)
+ read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
+ & qinfrag(i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
- & ifrag(2,i,iset), qinfrag(i,iset)
+ & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
+ & ifrag(2,i,iset1), qinfrag(i,iset1)
enddo
do i=1,npair
- read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
- & qinpair(i,iset)
+ read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
+ & qinpair(i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
- & ipair(2,i,iset), qinpair(i,iset)
+ & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
+ & ipair(2,i,iset1), qinpair(i,iset1)
enddo
do i=1,nfrag_back
- read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
- & wfrag_back(3,i,iset),
- & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
+ & wfrag_back(3,i,iset1),
+ & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
- & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ & write(iout,*) "A",i,wfrag_back(1,i,iset1),
+ & wfrag_back(2,i,iset1),
+ & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
+ & ifrag_back(2,i,iset1)
enddo
enddo
return
character*2 kic2
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
- integer ki, i, j, k, l
+ integer ki, i, j, k, l, ii_in_use(maxdim)
logical lprn /.true./
+ integer ilen
+ external ilen
c
c FP - Nov. 2014 Temporary specifications for new vars
c
c521 continue
c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+ ii_in_use(ii)=0
+ enddo
+ enddo
+
do k=1,constr_homology
read(inp,'(a)') pdbfile
c if(me.eq.king.or. .not. out1file)
c write (iout,'(2a)') 'PDB data will be read from file ',
c & pdbfile(:ilen(pdbfile))
+ 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)') 'Error opening PDB file.'
+ 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
+ & pdbfile(:ilen(pdbfile))
stop
34 continue
c print *,'Begin reading pdb data'
c enddo
c 1401 continue
c close (ientin)
- if (waga_dist.gt.0.0d0) then
+ if (waga_dist.ne.0.0d0) then
ii=0
do i = nnt,nct-2 ! right? without parallel.
do j=i+2,nct ! right?
c do j=i+2,nres ! ibid
c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
c do j=i+2,nct ! ibid
+
+ if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then
+
ii=ii+1
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+
c write (iout,*) "k",k
c write (iout,*) "i",i," j",j," constr_homology",
c & constr_homology
sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
else
+#ifdef OLDSIGMA
sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
& dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
+#else
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+#endif
c Following expr replaced by a positive exp argument
c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
c
c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
+ else
+ ii=ii+1
+ l_homo(k,ii)=.false.
+ endif
enddo
c read (ientin,*) sigma_odl(k,ii) ! 1st variant
enddo
c lim_odl=ii
c if (constr_homology.gt.0) call homology_partition
+ lim_odl=ii
endif
c
c Theta, dihedral and SC retraints
c Instead of res sim other local measure of b/b str reliability possible
sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
- if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
-c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
enddo
+ lim_dih=nct-nnt-2
endif
if (waga_theta.gt.0.0d0) then
c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
c rescore(k,i-2) ! right expression ?
c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
- if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
enddo
endif
+ lim_theta=nct-nnt-1
if (waga_d.gt.0.0d0) then
c open (ientin,file=tpl_k_sigma_d,status='old')
if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1 continue
enddo
+ lim_xx=nct-nnt+1
endif
close(ientin)
enddo
- if (waga_dist.gt.0.0d0) lim_odl=ii
+c
+c remove distance restraints not used in any model from the list
+c shift data in all arrays
+c
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i=nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+ if (ii_in_use(ii).eq.0) then
+ do ki=ii,lim_odl-1
+ ires_homo(ki)=ires_homo(ki+1)
+ jres_homo(ki)=jres_homo(ki+1)
+ ii_in_use(ki)=ii_in_use(ki+1)
+ do k=1,constr_homology
+ odl(k,ki)=odl(k,ki+1)
+ sigma_odl(k,ki)=sigma_odl(k,ki+1)
+ l_homo(k,ki)=l_homo(k,ki+1)
+ enddo
+ enddo
+ ii=ii-1
+ lim_odl=lim_odl-1
+ endif
+ enddo
+ enddo
+ endif
if (constr_homology.gt.0) call homology_partition
if (constr_homology.gt.0) call init_int_table
cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
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,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
- & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
+ write(iout,'(3i5,10(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,lim_dih
projects / unres.git / blobdiff