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,'NCUT',ncut,0)
+ call readi(controlcard,'NCLUST',nclust,5)
call readi(controlcard,'NSTART',nstart,0)
call readi(controlcard,'NEND',nend,0)
call reada(controlcard,'ECUT',ecut,10.0d0)
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)
+ if (ncut.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,'CONSTR_HOMOL',constr_homology,0)
write (iout,*) "with_homology_constr ",with_dihed_constr,
& " CONSTR_HOMOLOGY",constr_homology
+ print_homology_restraints=
+ & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
+ print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
+ print_homology_models=
+ & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
+#ifdef AIX
+ call flush_(iout)
+#else
call flush(iout)
+#endif
if (min_var) iopt=1
return
end
bad(i,2)=scalscp*bad(i,2)
enddo
+#ifdef AIX
+ call flush_(iout)
+#else
call flush(iout)
+#endif
print *,'indpdb=',indpdb,' pdbref=',pdbref
C Read sequence if not taken from the pdb file.
do i=1,nres
itype(i)=rescode(i,sequence(i),iscode)
enddo
+ if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
+ write (iout,*)
+ & "Glycine is the first full residue, initial dummy deleted"
+ do i=1,nres
+ itype(i)=itype(i+1)
+ enddo
+ nres=nres-1
+ endif
+ if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
+ write (iout,*)
+ & "Glycine is the last full residue, terminal dummy deleted"
+ nres=nres-1
+ endif
print *,nres
print '(20i4)',(itype(i),i=1,nres)
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
+#ifdef DEBUG
write (iout,*) "Calling init_dfa_vars"
- call flush(iout)
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+#endif
call init_dfa_vars
+#ifdef DEBUG
write (iout,*) 'init_dfa_vars finished!'
- call flush(iout)
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+#endif
call read_dfa_info
+#ifdef DEBUG
write (iout,*) 'read_dfa_info finished!'
- call flush(iout)
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+#endif
endif
if (constr_homology.gt.0) then
- call read_constr_homology
+ call read_constr_homology(print_homology_restraints)
endif
c if (pdbref) then
enddo
enddo
endif
- call contact(.true.,ncont_ref,icont_ref)
+ call contact(print_contact_map,ncont_ref,icont_ref)
endif
c Read distance restraints
if (constr_dist.gt.0) then
do i=1,npair_
write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
enddo
+#ifdef AIX
+ call flush_(iout)
+#else
call flush(iout)
+#endif
if (.not.refstr .and. nfrag_.gt.0) then
write (iout,*)
& "ERROR: no reference structure to compute distance restraints"
& " but distance restraints between pairs requested"
stop
endif
+#ifdef AIX
+ call flush_(iout)
+#else
call flush(iout)
+#endif
do i=1,nfrag_
if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
if (ifrag_(2,i).gt.nstart_sup+nsup-1)
& ifrag_(2,i)=nstart_sup+nsup-1
c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
- call flush(iout)
+c call flush(iout)
if (wfrag_(i).gt.0.0d0) then
do j=ifrag_(1,i),ifrag_(2,i)-1
do k=j+1,ifrag_(2,i)
write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
& i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
enddo
+#ifdef AIX
+ call flush_(iout)
+#else
call flush(iout)
+#endif
return
end
c====-------------------------------------------------------------------
- subroutine read_constr_homology
+ subroutine read_constr_homology(lprn)
include 'DIMENSIONS'
#ifdef MPI
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
integer ki, i, j, k, l
- logical lprn /.true./
+ logical lprn
+ logical unres_pdb
c
c FP - Nov. 2014 Temporary specifications for new vars
c
c
c
c Alternative: reading from input
+#ifdef DEBUG
write (iout,*) "BEGIN READ HOMOLOGY INFO"
+#ifdef AIX
+ call flush_(iout)
+#else
call flush(iout)
+#endif
+#endif
call card_concat(controlcard)
call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
call readi(controlcard,"HOMOL_NSET",homol_nset,1)
if (homol_nset.gt.1)then
+ call readi(controlcard,"ISET",iset,1)
call card_concat(controlcard)
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 "
-#ifdef DEBUG
- homol_nset=1
- call reada(controlcard,"WAGA_HOMOLOGY",waga_homology(1),1.0d0)
-#endif
- endif
- iset=mod(kolor,homol_nset)+1
else
- iset=1
- waga_homology(1)=1.0
+ iset=1
+ waga_homology(1)=1.0
endif
c
- write(iout,*) "read_constr_homology"
+#ifdef DEBUG
+ write(iout,*) "read_constr_homology iset",iset
write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
+#ifdef AIX
+ call flush_(iout)
+#else
call flush(iout)
-
-
+#endif
+#endif
cd write (iout,*) "nnt",nnt," nct",nct
cd call flush(iout)
do k=1,constr_homology
read(inp,'(a)') pdbfile
- write (iout,*) "k ",k," pdbfile ",pdbfile
+c write (iout,*) "k ",k," pdbfile ",pdbfile
c Next stament causes error upon compilation (?)
c if(me.eq.king.or. .not. out1file)
c write (iout,'(2a)') 'PDB data will be read from file ',
write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
& (crefjlee(j,i+nres),j=1,3)
enddo
-#endif
write (iout,*) "READ HOMOLOGY INFO"
write (iout,*) "read_constr_homology x: after reading pdb file"
write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
write (iout,*) "waga_theta",waga_theta
write (iout,*) "waga_d",waga_d
write (iout,*) "dist_cut",dist_cut
- call flush(iout)
+#endif
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
c
c Distance restraints
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
- sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
+#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 & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
c & "rescore(",k,i-2,") =",rescore(k,i-2),
c & "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) ! right expression ?
+ sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)+rescore(k,i-3))/4.0 ! right expression ?
c
c write (iout,*) "Raw sigmas for dihedral angle restraints"
c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
c & "rescore(",k,i-1,") =",rescore(k,i-1),
c & "rescore(",k,i-2,") =",rescore(k,i-2)
c read (ientin,*) sigma_theta(k,i) ! 1st variant
- sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
- & rescore(k,i-2) ! right expression ?
+ sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2))/3.0 ! right expression ?
sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*