time_Bcastw=time_Bcastw+MPI_Wtime()-time00
c call chainbuild_cart
endif
-c print *,'Processor',myrank,' calling etotal ipot=',ipot
+c write(iout,*) 'Processor',myrank,' calling etotal ipot=',ipot
c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
#else
c if (modecalc.eq.12.or.modecalc.eq.14) then
if (constr_homology.ge.1) then
call e_modeller(ehomology_constr)
+ print *,'iset=',iset,'me=',me,ehomology_constr,
+ & 'Processor',fg_rank,' CG group',kolor,
+ & ' absolute rank',MyRank
else
ehomology_constr=0.0d0
endif
C If performing constraint dynamics, call the constraint energy
C after the equilibration time
if(usampl.and.totT.gt.eq_time) then
+c write (iout,*) "CALL TO ECONSTR_BACK"
call EconstrQ
call Econstr_back
else
common /sccalc/ time11,time12,time112,theti,it,nlobit
delta=0.02d0*pi
escloc=0.0D0
+c write(iout,*) "ESC: loc_start",loc_start," loc_end",loc_end
do i=loc_start,loc_end
costtab(i+1) =dcos(theta(i+1))
sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
integer nnn, i, j, k, ki, irec, l
integer katy, odleglosci, test7
real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template)
+ real*8 Eval,Erot
real*8 distance(max_template),distancek(max_template),
& min_odl,godl(max_template),dih_diff(max_template)
+c
+c FP - 30/10/2014 Temporary specifications for homology restraints
+c
+ double precision utheta_i,gutheta_i,sum_gtheta,sum_sgtheta,
+ & sgtheta
+ double precision, dimension (maxres) :: guscdiff,usc_diff
+ double precision, dimension (max_template) ::
+ & gtheta,dscdiff,uscdiffk,guscdiff2,guscdiff3,
+ & theta_diff
+c
+
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.GEO'
include 'COMMON.IOUNITS'
include 'COMMON.MD'
include 'COMMON.CONTROL'
+c
+c From subroutine Econstr_back
+c
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+c
do i=1,19
c Pseudo-energy and gradient from homology restraints (MODELLER-like
c function)
C AL 5/2/14 - Introduce list of restraints
+c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs start -------"
+#endif
do ii = link_start_homo,link_end_homo
i = ires_homo(ii)
j = jres_homo(ii)
dij=dist(i,j)
+c write (iout,*) "dij(",i,j,") =",dij
do k=1,constr_homology
distance(k)=odl(k,ii)-dij
- distancek(k)=
- & 0.5d0*waga_dist(iset)*distance(k)**2*sigma_odl(k,ii)
+c write (iout,*) "distance(",k,") =",distance(k)
+ distancek(k)=0.5d0*distance(k)**2*sigma_odl(k,ii) ! waga_dist rmvd from Gaussian argument
+c write (iout,*) "sigma_odl(",k,ii,") =",sigma_odl(k,ii)
+c write (iout,*) "distancek(",k,") =",distancek(k)
+c distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii)
enddo
min_odl=minval(distancek)
+c write (iout,* )"min_odl",min_odl
#ifdef DEBUG
write (iout,*) "ij dij",i,j,dij
write (iout,*) "distance",(distance(k),k=1,constr_homology)
odleg2=0.0d0
do k=1,constr_homology
c Nie wiem po co to liczycie jeszcze raz!
- c odleg3=-waga_dist*((distance(i,j,k)**2)/
+ c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
c & (2*(sigma_odl(i,j,k))**2))
godl(k)=dexp(-distancek(k)+min_odl)
odleg2=odleg2+godl(k)
ccc & "sigma_odl(i,j,k)=", sigma_odl(i,j,k)
enddo
+c write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+c write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
#ifdef DEBUG
- write (iout,*) "godl",(godl(k),k=1,constr_homology)
- write (iout,*) "ii i j",ii,i,j," odleg2",odleg2
+ write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+ write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
#endif
odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+c write (iout,*) "odleg",odleg ! sum of -ln-s
c Gradient
sum_godl=odleg2
- sum_sgodl=0.0
+ sum_sgodl=0.0d0
do k=1,constr_homology
c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
-c & *waga_dist(iset))+min_odl
- sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist(iset)
+c & *waga_dist)+min_odl
+c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
+ sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd
sum_sgodl=sum_sgodl+sgodl
c sgodl2=sgodl2+sgodl
c write(iout,*) i, j, k, "TEST K"
enddo
- grad_odl3=sum_sgodl/(sum_godl*dij)
++ if (homol_nset.gt.1)then
++ grad_odl3=waga_dist1(iset)*sum_sgodl/(sum_godl*dij)
++ else
+ grad_odl3=waga_dist*sum_sgodl/(sum_godl*dij)
++ endif
+c grad_odl3=sum_sgodl/(sum_godl*dij)
c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2"
ghpbc(jik,j)=ghpbc(jik,j)-ggodl
ccc write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl,
ccc & ghpbc(jik,i+1), ghpbc(jik,j+1)
-
+c if (i.eq.25.and.j.eq.27) then
+c write(iout,*) "jik",jik,"i",i,"j",j
+c write(iout,*) "sum_sgodl",sum_sgodl,"sgodl",sgodl
+c write(iout,*) "grad_odl3",grad_odl3
+c write(iout,*) "c(",jik,i,")",c(jik,i),"c(",jik,j,")",c(jik,j)
+c write(iout,*) "ggodl",ggodl
+c write(iout,*) "ghpbc(",jik,i,")",
+c & ghpbc(jik,i),"ghpbc(",jik,j,")",
+c & ghpbc(jik,j)
+c endif
enddo
ccc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=",
ccc & dLOG(odleg2),"-odleg=", -odleg
- enddo ! ii
+ enddo ! ii-loop for dist
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs end -------"
+c if (waga_angle.eq.1.0d0 .or. waga_theta.eq.1.0d0 .or.
+c & waga_d.eq.1.0d0) call sum_gradient
+#endif
c Pseudo-energy and gradient from dihedral-angle restraints from
c homology templates
c write (iout,*) "End of distance loop"
c call flush(iout)
kat=0.0d0
c write (iout,*) idihconstr_start_homo,idihconstr_end_homo
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs start -------"
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ write (iout,*) "gloc_init(",i,icg,")",gloc(i,icg)
+ enddo
+#endif
do i=idihconstr_start_homo,idihconstr_end_homo
kat2=0.0d0
c betai=beta(i,i+1,i+2,i+3)
betai = phi(i+3)
+c write (iout,*) "betai =",betai
do k=1,constr_homology
dih_diff(k)=pinorm(dih(k,i)-betai)
+c write (iout,*) "dih_diff(",k,") =",dih_diff(k)
c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
c & -(6.28318-dih_diff(i,k))
c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
c & 6.28318+dih_diff(i,k)
- kat3=-0.5d0*waga_angle(iset)*dih_diff(k)**2*sigma_dih(k,i)
+ kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
gdih(k)=dexp(kat3)
kat2=kat2+gdih(k)
c write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
c write(*,*)""
enddo
+c write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps
+c write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps
#ifdef DEBUG
write (iout,*) "i",i," betai",betai," kat2",kat2
write (iout,*) "gdih",(gdih(k),k=1,constr_homology)
#endif
if (kat2.le.1.0d-14) cycle
kat=kat-dLOG(kat2/constr_homology)
+c write (iout,*) "kat",kat ! sum of -ln-s
ccc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=",
ccc & dLOG(kat2), "-kat=", -kat
c ----------------------------------------------------------------------
sum_gdih=kat2
- sum_sgdih=0.0
+ sum_sgdih=0.0d0
do k=1,constr_homology
- sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle(iset)
+ sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd
+c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
sum_sgdih=sum_sgdih+sgdih
enddo
- grad_dih3=sum_sgdih/sum_gdih
+c grad_dih3=sum_sgdih/sum_gdih
- grad_dih3=waga_angle*sum_sgdih/sum_gdih
++ if (homol_nset.gt.1)then
++ grad_dih3=waga_angle1(iset)*sum_sgdih/sum_gdih
++ else
++ grad_dih3=waga_angle*sum_sgdih/sum_gdih
++ endif
c write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
ccc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
ccc & gloc(nphi+i-3,icg)
gloc(i,icg)=gloc(i,icg)+grad_dih3
+c if (i.eq.25) then
+c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
+c endif
ccc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3,
ccc & gloc(nphi+i-3,icg)
+ enddo ! i-loop for dih
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs end -------"
+#endif
+
+c Pseudo-energy and gradient for theta angle restraints from
+c homology templates
+c FP 01/15 - inserted from econstr_local_test.F, loop structure
+c adapted
+
+c
+c For constr_homology reference structures (FP)
+c
+c Uconst_back_tot=0.0d0
+ Eval=0.0d0
+ Erot=0.0d0
+c Econstr_back legacy
+ do i=1,nres
+c do i=ithet_start,ithet_end
+ dutheta(i)=0.0d0
+c enddo
+c do i=loc_start,loc_end
+ do j=1,3
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
+ enddo
enddo
+c
+c do iref=1,nref
+c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c write (iout,*) "waga_theta",waga_theta
+ if (waga_theta.gt.0.0d0) then
+#ifdef DEBUG
+ write (iout,*) "usampl",usampl
+ write(iout,*) "------- theta restrs start -------"
+c do i=ithet_start,ithet_end
+c write (iout,*) "gloc_init(",nphi+i,icg,")",gloc(nphi+i,icg)
+c enddo
+#endif
+c write (iout,*) "maxres",maxres,"nres",nres
+
+ do i=ithet_start,ithet_end
+c
+c do i=1,nfrag_back
+c ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset)
+c
+c Deviation of theta angles wrt constr_homology ref structures
+c
+ utheta_i=0.0d0 ! argument of Gaussian for single k
+ gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+c do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop
+c over residues in a fragment
+c write (iout,*) "theta(",i,")=",theta(i)
+ do k=1,constr_homology
+c
+c dtheta_i=theta(j)-thetaref(j,iref)
+c dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing
+ theta_diff(k)=thetatpl(k,i)-theta(i)
+c
+ utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
+c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
+ gtheta(k)=dexp(utheta_i) ! + min_utheta_i?
+ gutheta_i=gutheta_i+dexp(utheta_i) ! Sum of Gaussians (pk)
+c Gradient for single Gaussian restraint in subr Econstr_back
+c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
+c
+ enddo
+c write (iout,*) "gtheta",(gtheta(k),k=1,constr_homology) ! exps
+c write (iout,*) "i",i," gutheta_i",gutheta_i ! sum of exps
+c
+c Gradient for multiple Gaussian restraint
+ sum_gtheta=gutheta_i
+ sum_sgtheta=0.0d0
+ do k=1,constr_homology
+c New generalized expr for multiple Gaussian from Econstr_back
+ sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i) ! waga_theta rmvd
+c
+c sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form?
+ sum_sgtheta=sum_sgtheta+sgtheta ! cum variable
+ enddo
+c grad_theta3=sum_sgtheta/sum_gtheta 1/*theta(i)? s. line below
+c grad_theta3=sum_sgtheta/sum_gtheta
+c
+c Final value of gradient using same var as in Econstr_back
+ dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta
+c dutheta(i)=sum_sgtheta/sum_gtheta
+c
+c Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight
+ Eval=Eval-dLOG(gutheta_i/constr_homology)
+c write (iout,*) "utheta(",i,")=",utheta(i) ! -ln of sum of exps
+c write (iout,*) "Uconst_back",Uconst_back ! sum of -ln-s
+c Uconst_back=Uconst_back+utheta(i)
+ enddo ! (i-loop for theta)
+#ifdef DEBUG
+ write(iout,*) "------- theta restrs end -------"
+#endif
+ endif
+c
+c Deviation of local SC geometry
+c
+c Separation of two i-loops (instructed by AL - 11/3/2014)
+c
+c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c write (iout,*) "waga_d",waga_d
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs start -------"
+ write (iout,*) "Initial duscdiff,duscdiffx"
+ do i=loc_start,loc_end
+ write (iout,*) i,(duscdiff(jik,i),jik=1,3),
+ & (duscdiffx(jik,i),jik=1,3)
+ enddo
+#endif
+ do i=loc_start,loc_end
+ usc_diff_i=0.0d0 ! argument of Gaussian for single k
+ guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+c do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy
+c write(iout,*) "xxtab, yytab, zztab"
+c write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i)
+ do k=1,constr_homology
+c
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+c Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+c write(iout,*) "dxx, dyy, dzz"
+c write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+c
+ usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument
+c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
+c uscdiffk(k)=usc_diff(i)
+ guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
+ guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk)
+c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
+c & xxref(j),yyref(j),zzref(j)
+ enddo
+c
+c Gradient
+c
+c Generalized expression for multiple Gaussian acc to that for a single
+c Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014)
+c
+c Original implementation
+c sum_guscdiff=guscdiff(i)
+c
+c sum_sguscdiff=0.0d0
+c do k=1,constr_homology
+c sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d?
+c sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff
+c sum_sguscdiff=sum_sguscdiff+sguscdiff
+c enddo
+c
+c Implementation of new expressions for gradient (Jan. 2015)
+c
+c grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !?
+ do k=1,constr_homology
+c
+c New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong
+c before. Now the drivatives should be correct
+c
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+c Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+c
+c New implementation
+c
+ sum_guscdiff=guscdiff2(k)*!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong!
+ & sigma_d(k,i) ! for the grad wrt r'
+c sum_sguscdiff=sum_sguscdiff+sum_guscdiff
+c
+c
+c New implementation
+ sum_guscdiff = waga_d*sum_guscdiff
+ do jik=1,3
+ duscdiff(jik,i-1)=duscdiff(jik,i-1)+
+ & sum_guscdiff*(dXX_C1tab(jik,i)*dxx+
+ & dYY_C1tab(jik,i)*dyy+dZZ_C1tab(jik,i)*dzz)/guscdiff(i)
+ duscdiff(jik,i)=duscdiff(jik,i)+
+ & sum_guscdiff*(dXX_Ctab(jik,i)*dxx+
+ & dYY_Ctab(jik,i)*dyy+dZZ_Ctab(jik,i)*dzz)/guscdiff(i)
+ duscdiffx(jik,i)=duscdiffx(jik,i)+
+ & sum_guscdiff*(dXX_XYZtab(jik,i)*dxx+
+ & dYY_XYZtab(jik,i)*dyy+dZZ_XYZtab(jik,i)*dzz)/guscdiff(i)
+c
+#ifdef DEBUG
+ write(iout,*) "jik",jik,"i",i
+ write(iout,*) "dxx, dyy, dzz"
+ write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+ write(iout,*) "guscdiff2(",k,")",guscdiff2(k)
+c write(iout,*) "sum_sguscdiff",sum_sguscdiff
+cc write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i)
+c write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i)
+c write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i)
+c write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i)
+c write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i)
+c write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i)
+c write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i)
+c write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i)
+c write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i)
+c write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1)
+c write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i)
+c write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i)
+c endif
+#endif
+ enddo
+ enddo
+c
+c uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required?
+c usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ?
+c
+c write (iout,*) i," uscdiff",uscdiff(i)
+c
+c Put together deviations from local geometry
+
+c Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+
+c & wfrag_back(3,i,iset)*uscdiff(i)
+ Erot=Erot-dLOG(guscdiff(i)/constr_homology)
+c write (iout,*) "usc_diff(",i,")=",usc_diff(i) ! -ln of sum of exps
+c write (iout,*) "Uconst_back",Uconst_back ! cum sum of -ln-s
+c Uconst_back=Uconst_back+usc_diff(i)
+c
+c Gradient of multiple Gaussian restraint (FP - 04/11/2014 - right?)
+c
+c New implment: multiplied by sum_sguscdiff
+c
+
+ enddo ! (i-loop for dscdiff)
+
+c endif
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs end -------"
+ write (iout,*) "------ After SC loop in e_modeller ------"
+ do i=loc_start,loc_end
+ write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
+ write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
+ enddo
+ if (waga_theta.eq.1.0d0) then
+ write (iout,*) "in e_modeller after SC restr end: dutheta"
+ do i=ithet_start,ithet_end
+ write (iout,*) i,dutheta(i)
+ enddo
+ endif
+ if (waga_d.eq.1.0d0) then
+ write (iout,*) "e_modeller after SC loop: duscdiff/x"
+ do i=1,nres
+ write (iout,*) i,(duscdiff(j,i),j=1,3)
+ write (iout,*) i,(duscdiffx(j,i),j=1,3)
+ enddo
+ endif
+#endif
c Total energy from homology restraints
#ifdef DEBUG
write (iout,*) "odleg",odleg," kat",kat
#endif
- ehomology_constr=odleg+kat
+c
+c Addition of energy of theta angle and SC local geom over constr_homologs ref strs
+c
+c ehomology_constr=odleg+kat
- ehomology_constr=waga_dist*odleg+waga_angle*kat+waga_theta*Eval
++ if (homol_nset.gt.1)then
++ ehomology_constr=waga_dist1(iset)*odleg+waga_angle1(iset)*kat+waga_theta*Eval
++ & +waga_d*Erot
++ else
++ ehomology_constr=waga_dist*odleg+waga_angle*kat+waga_theta*Eval
+ & +waga_d*Erot
++ endif
+c write (iout,*) "odleg",odleg," kat",kat," Uconst_back",Uconst_back
+c write (iout,*) "ehomology_constr",ehomology_constr
+c ehomology_constr=odleg+kat+Uconst_back
return
-
+c
+c FP 01/15 end
+c
748 format(a8,f12.3,a6,f12.3,a7,f12.3)
747 format(a12,i4,i4,i4,f8.3,f8.3)
746 format(a12,i4,i4,i4,f8.3,f8.3,f8.3)
endif
call reada(controlcard,"Q_NP",Q_np,0.1d0)
usampl = index(controlcard,"USAMPL").gt.0
-
mdpdb = index(controlcard,"MDPDB").gt.0
call reada(controlcard,"T_BATH",t_bath,300.0d0)
call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
if (constr_homology.gt.0) then
call read_constr_homology
+ else
+ homol_nset=0
endif
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.MD'
+ include 'COMMON.CONTROL'
open(irest2,file=rest2name,status='unknown')
read(irest2,*) totT,EK,potE,totE,t_bath
do i=1,2*nres
do i=1,2*nres
read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
- if(usampl) then
+ if(usampl.or.homol_nset.gt.1) then
read (irest2,*) iset
endif
close(irest2)
include 'COMMON.MD'
include 'COMMON.GEO'
include 'COMMON.INTERACT'
- double precision odl_temp,sigma_odl_temp
- common /przechowalnia/ odl_temp(maxres,maxres,max_template),
- & sigma_odl_temp(maxres,maxres,max_template)
+c
+c For new homol impl
+c
+ include 'COMMON.VAR'
+c
+
+c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
+c & dist_cut
+c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+c & sigma_odl_temp(maxres,maxres,max_template)
character*2 kic2
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
- character*3200 controlcard1
integer ki, i, j, k, l
logical lprn /.true./
-
+c
+c FP - Nov. 2014 Temporary specifications for new vars
+c
+ double precision rescore_tmp,x12,y12,z12
+ double precision, dimension (max_template,maxres) :: rescore
+ character*24 pdbfile,tpl_k_rescore
+c -----------------------------------------------------------------
+c Reading multiple PDB ref structures and calculation of retraints
+c not using pre-computed ones stored in files model_ki_{dist,angle}
+c FP (Nov., 2014)
+c -----------------------------------------------------------------
+c
+c
+c Alternative: reading from input
call card_concat(controlcard)
+ 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 readi(controlcard,"HOMOL_NSET",homol_nset,1)
+ if (homol_nset.gt.1)then
+ call card_concat(controlcard)
- read(controlcard,*) (waga_dist(i),i=1,homol_nset)
++ read(controlcard,*) (waga_dist1(i),i=1,homol_nset)
+ call card_concat(controlcard)
- read(controlcard,*) (waga_angle(i),i=1,homol_nset)
++ read(controlcard,*) (waga_angle1(i),i=1,homol_nset)
+ write(iout,*) "iset distance_weight angle_weight"
+ do i=1,homol_nset
- write(iout,*) i,waga_dist(i),waga_angle(i)
++ write(iout,*) i,waga_dist1(i),waga_angle1(i)
+ enddo
+ iset=mod(kolor,homol_nset)+1
+ else
+ iset=1
- call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
- call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
++c call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
++c call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
+ endif
+
+ write (iout,*) "nnt",nnt," nct",nct
+ call flush(iout)
++
++
lim_odl=0
lim_dih=0
- do i=1,nres
- do j=i+2,nres
- do ki=1,constr_homology
- sigma_odl_temp(i,j,ki)=0.0d0
- odl_temp(i,j,ki)=0.0d0
+c
+c New
+c
+ lim_theta=0
+ lim_xx=0
+c
+c Reading HM global scores (prob not required)
+c
+c open (4,file="HMscore")
+c do k=1,constr_homology
+c read (4,*,end=521) hmscore_tmp
+c hmscore(k)=hmscore_tmp ! Another transformation can be used
+c write(*,*) "Model", k, ":", hmscore(k)
+c enddo
+c521 continue
+
+c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+ do k=1,constr_homology
+
+ read(inp,'(a)') 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 ',
+c & pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a)') 'Error opening PDB file.'
+ stop
+ 34 continue
+c print *,'Begin reading pdb data'
+c
+c Files containing res sim or local scores (former containing sigmas)
+c
+
+ write(kic2,'(bz,i2.2)') k
+
+ tpl_k_rescore="template"//kic2//".sco"
+c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
+c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
+c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
+c tpl_k_sigma_d="template"//kic2//".sigma_d"
+
+ call readpdb
+c
+c Distance restraints
+c
+c ... --> odl(k,ii)
+C Copy the coordinates from reference coordinates (?)
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=cref(j,i)
+c write (iout,*) "c(",j,i,") =",c(j,i)
enddo
enddo
- enddo
- do i=1,nres-3
- do ki=1,constr_homology
- dih(ki,i)=0.0d0
- sigma_dih(ki,i)=0.0d0
- enddo
- enddo
- do ki=1,constr_homology
- write(kic2,'(i2)') ki
- if (ki.le.9) kic2="0"//kic2(2:2)
+c
+c From read_dist_constr (commented out 25/11/2014 <-> res sim)
+c
+c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
+ open (ientin,file=tpl_k_rescore,status='old')
+ do irec=1,maxdim ! loop for reading res sim
+ if (irec.eq.1) then
+ rescore(k,irec)=0.0d0
+ goto 1301
+ endif
+ read (ientin,*,end=1401) rescore_tmp
+c 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
+c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
+ 1301 continue
+ enddo
+ 1401 continue
+ close (ientin)
+c open (ientin,file=tpl_k_sigma_odl,status='old')
+c do irec=1,maxdim ! loop for reading sigma_odl
+c read (ientin,*,end=1401) i, j,
+c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
+c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
+c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
+c enddo
+c 1401 continue
+c close (ientin)
+ if (waga_dist.gt.0.0d0) then
+ ii=0
+ do i = nnt,nct-2 ! right? without parallel.
+ do j=i+2,nct ! right?
+c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
+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
+ ii=ii+1
+c write (iout,*) "k",k
+c write (iout,*) "i",i," j",j," constr_homology",
+c & constr_homology
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+c
+c Attempt to replace dist(i,j) by its definition in ...
+c
+ 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)
+ odl(k,ii)=distal
+c
+c odl(k,ii)=dist(i,j)
+c write (iout,*) "dist(",i,j,") =",dist(i,j)
+c write (iout,*) "distal = ",distal
+c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,j,") =",rescore(k,j)
+c
+c Calculation of sigma from res sim
+c
+c if (odl(k,ii).le.6.0d0) then
+c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
+c Other functional forms possible depending on odl(k,ii), eg.
+c
+ if (odl(k,ii).le.dist_cut) then
+ 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
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
- model_ki_dist="model"//kic2//".dist"
- model_ki_angle="model"//kic2//".angle"
- open (ientin,file=model_ki_dist,status='old')
- do irec=1,maxdim !petla do czytania wiezow na odleglosc
- read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
- & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
- odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
- sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
- & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
- enddo
- 1401 continue
- close (ientin)
- open (ientin,file=model_ki_angle,status='old')
- do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
- read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
- & sigma_dih(ki,i+nnt-1)
- if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
- sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
- enddo
- 1402 continue
- close (ientin)
- enddo
- ii=0
- write (iout,*) "nnt",nnt," nct",nct
- do i=nnt,nct-2
- do j=i+2,nct
- ki=1
-c write (iout,*) "i",i," j",j," constr_homology",constr_homology
- do while (ki.le.constr_homology .and.
- & sigma_odl_temp(i,j,ki).le.0.0d0)
-c write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
- ki=ki+1
+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 sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
+c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+ endif
+c
+ sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
+c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
+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
+ enddo
+c read (ientin,*) sigma_odl(k,ii) ! 1st variant
enddo
-c write (iout,*) "ki",ki
- if (ki.gt.constr_homology) cycle
- ii=ii+1
- ires_homo(ii)=i
- jres_homo(ii)=j
- do ki=1,constr_homology
- odl(ki,ii)=odl_temp(i,j,ki)
- sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
- enddo
- enddo
+c lim_odl=ii
+c if (constr_homology.gt.0) call homology_partition
+ endif
+c
+c Theta, dihedral and SC retraints
+c
+ if (waga_angle.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_dih,status='old')
+c do irec=1,maxres-3 ! loop for reading sigma_dih
+c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
+c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
+c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_dih(k,i+nnt-1)
+c enddo
+c1402 continue
+c close (ientin)
+ do i = nnt+3,nct ! right? without parallel.
+c do i=1,nres ! alternative for bounds acc to readpdb?
+c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
+c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
+ dih(k,i)=phiref(i) ! right?
+c read (ientin,*) sigma_dih(k,i) ! original variant
+c write (iout,*) "dih(",k,i,") =",dih(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,i-1,") =",rescore(k,i-1),
+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 ?
+c
+c write (iout,*) "Raw sigmas for dihedral angle restraints"
+c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
+c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
+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
+ endif
+
+ if (waga_theta.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_theta,status='old')
+c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
+c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
+c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_theta(k,i+nnt-1)
+c enddo
+c1403 continue
+c close (ientin)
+
+ do i = nnt+2,nct ! right? without parallel.
+c do i = i=1,nres ! alternative for bounds acc to readpdb?
+c do i=ithet_start,ithet_end ! with FG parallel.
+ thetatpl(k,i)=thetaref(i)
+c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(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)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+
+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
+
+ if (waga_d.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_d,status='old')
+c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
+c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
+c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_d(k,i+nnt-1)
+c enddo
+c1404 continue
+ close (ientin)
+
+ do i = nnt,nct ! right? without parallel.
+c do i=2,nres-1 ! alternative for bounds acc to readpdb?
+c do i=loc_start,loc_end ! with FG parallel.
+ if (itype(i).eq.10) goto 1 ! right?
+ xxtpl(k,i)=xxref(i)
+ yytpl(k,i)=yyref(i)
+ zztpl(k,i)=zzref(i)
+c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
+c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
+c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
+ sigma_d(k,i)=rescore(k,i) ! right expression ?
+ sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+
+c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
+c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
+c read (ientin,*) sigma_d(k,i) ! 1st variant
+ if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
+ 1 continue
+ enddo
+ endif
+ close(ientin)
enddo
- lim_odl=ii
+ if (waga_dist.gt.0.0d0) lim_odl=ii
if (constr_homology.gt.0) call homology_partition
+ if (constr_homology.gt.0) call init_int_table
+ write (iout,*) "homology_partition: lim_theta= ",lim_theta,
+ & "lim_xx=",lim_xx
+c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c
c Print restraints
+c
if (.not.lprn) return
+ write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
write (iout,*) "Distance restraints from templates"
do ii=1,lim_odl
- write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
+ 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)
enddo
write (iout,*) "Dihedral angle restraints from templates"
- do i=nnt,lim_dih
+ do i=nnt+3,lim_dih
write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
& rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
enddo
-c write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
-c write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
-
+ write (iout,*) "Virtual-bond angle restraints from templates"
+ do i=nnt+2,lim_theta
+ write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+ & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "SC restraints from templates"
+ do i=nnt,lim_xx
+ write(iout,'(i5,10(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
+c -----------------------------------------------------------------
return
end
c----------------------------------------------------------------------
projects / unres.git / commitdiff