From: Felipe Pineda Date: Thu, 12 Feb 2015 13:54:55 +0000 (+0100) Subject: 1st running version of UNRES HM by FP and AL X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=73ca05fce4b785367bd9600c24a7f026feff14bf;p=unres.git 1st running version of UNRES HM by FP and AL --- diff --git a/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe b/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe index 0929cef..16d090c 100755 Binary files a/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe and b/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe differ diff --git a/bin/unres/MD/unres_ifort_MPICH_GAB.exe b/bin/unres/MD/unres_ifort_MPICH_GAB.exe index bce0500..e523a19 100755 Binary files a/bin/unres/MD/unres_ifort_MPICH_GAB.exe and b/bin/unres/MD/unres_ifort_MPICH_GAB.exe differ diff --git a/source/unres/src_MD/COMMON.CONTROL b/source/unres/src_MD/COMMON.CONTROL index 9fce3c5..ce2c00d 100644 --- a/source/unres/src_MD/COMMON.CONTROL +++ b/source/unres/src_MD/COMMON.CONTROL @@ -1,6 +1,6 @@ integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad, & inprint,i2ndstr,mucadyn,constr_dist,constr_homology - real*8 waga_dist, waga_angle + real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec, & sideadd,lsecondary,read_cart,unres_pdb, & vdisulf,searchsc,lmuca,dccart,extconf,out1file, @@ -10,6 +10,6 @@ & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file, & constr_dist,gnorm_check,gradout,split_ene,constr_homology, - & waga_dist, waga_angle + & waga_dist, waga_angle, waga_theta, waga_d, dist_cut C... minim = .true. means DO minimization. C... energy_dec = .true. means print energy decomposition matrix diff --git a/source/unres/src_MD/COMMON.MD b/source/unres/src_MD/COMMON.MD index bd38d1b..b66ea2c 100644 --- a/source/unres/src_MD/COMMON.MD +++ b/source/unres/src_MD/COMMON.MD @@ -16,6 +16,14 @@ real*8 odl(max_template,maxdim),sigma_odl(max_template,maxdim), & dih(max_template,maxres),sigma_dih(max_template,maxres) +c +c Specification of new variables used in subroutine e_modeller +c modified by FP (Nov.,2014) + real*8 xxtpl(max_template,maxres),yytpl(max_template,maxres), + & zztpl(max_template,maxres),thetatpl(max_template,maxres), + & sigma_theta(max_template,maxres), + & sigma_d(max_template,maxres) +c integer ires_homo(maxdim),jres_homo(maxdim) @@ -36,6 +44,11 @@ & ifrag_back(3,maxfrag_back,maxprocs/20),ntime_split,ntime_split0, & maxtime_split,lim_odl,lim_dih,link_start_homo,link_end_homo, & idihconstr_start_homo,idihconstr_end_homo +c +c FP (30/10/2014) +c +c integer ithetaconstr_start_homo,ithetaconstr_end_homo +c integer nresn,nyosh,nnos double precision glogs,qmass,vlogs,xlogs logical large,print_compon,tbf,rest,reset_moment,reset_vel, @@ -47,7 +60,12 @@ & wfrag_back,nfrag_back,ifrag_back common /homrestr/ odl,dih,sigma_dih,sigma_odl, & lim_odl,lim_dih,ires_homo,jres_homo,link_start_homo, - & link_end_homo,idihconstr_start_homo,idihconstr_end_homo + & link_end_homo,idihconstr_start_homo,idihconstr_end_homo, +c +c FP (30/10/2014) +c + & xxtpl,yytpl,zztpl,thetatpl,sigma_theta,sigma_d +c common /qmeas/ qfrag,qpair,qinfrag,qinpair,wfrag,wpair,eq_time, & Ucdfrag,Ucdpair,dUdconst,dUdxconst,dqwol,dxqwol,Uconst, & iset,mset,nset,usampl,ifrag,ipair,npair,nfrag diff --git a/source/unres/src_MD/Makefile_MPICH_ifort b/source/unres/src_MD/Makefile_MPICH_ifort index d763388..ebdb428 100644 --- a/source/unres/src_MD/Makefile_MPICH_ifort +++ b/source/unres/src_MD/Makefile_MPICH_ifort @@ -10,6 +10,8 @@ FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include FFLAGS1 = -c -g -CA -CB -I$(INSTALL_DIR)/include FFLAGS2 = -c -g -O0 -I$(INSTALL_DIR)/include FFLAGSE = -c -O3 -ipo -opt_report -I$(INSTALL_DIR)/include +# FFLAGS = ${FFLAGS1} +# FFLAGSE = ${FFLAGS1} LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a diff --git a/source/unres/src_MD/checkder_p.F b/source/unres/src_MD/checkder_p.F index 4d0379e..7cf3ed6 100644 --- a/source/unres/src_MD/checkder_p.F +++ b/source/unres/src_MD/checkder_p.F @@ -395,6 +395,7 @@ c call int_from_cart1(.false.) if (.not.split_ene) then call etotal(energia1(0)) etot1=energia1(0) +c write (iout,*) "i",i," etot",etot," etot1",etot1 else !- split gradient call etotal_long(energia1(0)) @@ -441,6 +442,7 @@ c write (iout,*) if (.not.split_ene) then call etotal(energia1(0)) etot1=energia1(0) +c write (iout,*) "i",i," etot",etot," etot1",etot1 else !- split gradient call etotal_long(energia1(0)) @@ -686,7 +688,7 @@ cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar) call etotal(energia2(0)) etot2=energia2(0) gg(i)=(etot2-etot1)/aincr - write (iout,*) i,etot1,etot2 +c write (iout,*) i,etot1,etot2 x(i)=xi enddo write (iout,'(/2a)')' Variable Numerical Analytical', diff --git a/source/unres/src_MD/cinfo.f b/source/unres/src_MD/cinfo.f index b42001a..7cda99a 100644 --- a/source/unres/src_MD/cinfo.f +++ b/source/unres/src_MD/cinfo.f @@ -1,32 +1,31 @@ C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C -C 3 2 169 +C 3 2 304 subroutine cinfo include 'COMMON.IOUNITS' write(iout,*)'++++ Compile info ++++' - write(iout,*)'Version 3.2 build 169' - write(iout,*)'compiled Mon Jul 21 16:29:49 2014' - write(iout,*)'compiled by jal47@matrix.chem.cornell.edu' + write(iout,*)'Version 3.2 build 304' + write(iout,*)'compiled Thu Feb 5 13:01:33 2015' + write(iout,*)'compiled by felipe@piasek4' write(iout,*)'OS name: Linux ' - write(iout,*)'OS release: 2.6.34.9-69.fc13.x86_64 ' + write(iout,*)'OS release: 3.2.0-70-generic ' write(iout,*)'OS version:', - & ' #1 SMP Tue May 3 09:23:03 UTC 2011 ' + & ' #105-Ubuntu SMP Wed Sep 24 19:49:16 UTC 2014 ' write(iout,*)'flags:' write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...' - write(iout,*)'FC = ifort' - write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include' - write(iout,*)'FFLAGS1 = -c -w -g -d2 -CA -CB -I$(INSTALL_DIR)...' - write(iout,*)'FFLAGS2 = -c -w -g -O0 -I$(INSTALL_DIR)/include' - write(iout,*)'FFLAGS3 = -c -w -O3 -mp' - write(iout,*)'FFLAGSE = -c -w -O3 -ipo -ipo_obj -opt_report ...' - write(iout,*)'CC = cc' - write(iout,*)'CFLAGS = -DLINUX -DPGI -c' - write(iout,*)'OPT = -O3 -ip -w' + write(iout,*)'FC= ifort' + write(iout,*)'OPT = -O3 -ip ' + write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include ' + write(iout,*)'FFLAGS1 = -c -g -CA -CB -I$(INSTALL_DIR)/inclu...' + write(iout,*)'FFLAGS2 = -c -g -O0 -I$(INSTALL_DIR)/include ' + write(iout,*)'FFLAGSE = -c -O3 -ipo -opt_report -I$(INSTALL...' + write(iout,*)'FFLAGS = ${FFLAGS1}' + write(iout,*)'FFLAGSE = ${FFLAGS1}' write(iout,*)'LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdr...' write(iout,*)'ARCH = LINUX' write(iout,*)'PP = /lib/cpp -P' write(iout,*)'object = unres.o arcos.o cartprint.o chainbuild...' write(iout,*)'GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES ...' - write(iout,*)'GAB: BIN = ../bin/unres_ifort_MPICH-restr-DFA_G...' + write(iout,*)'GAB: BIN = ../../../bin/unres/MD/unres_ifort_MP...' write(iout,*)'E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNR...' write(iout,*)'E0LL2Y: BIN = ../../../bin/unres/MD/unres_ifort...' write(iout,*)'++++ End of compile info ++++' diff --git a/source/unres/src_MD/econstr_local.F b/source/unres/src_MD/econstr_local.F index f11acfb..88de310 100644 --- a/source/unres/src_MD/econstr_local.F +++ b/source/unres/src_MD/econstr_local.F @@ -19,14 +19,6 @@ c MD with umbrella_sampling using Wolyne's distance measure as a constraint include 'COMMON.NAMES' include 'COMMON.TIME1' Uconst_back=0.0d0 - do i=1,nres - dutheta(i)=0.0d0 - dugamma(i)=0.0d0 - do j=1,3 - duscdiff(j,i)=0.0d0 - duscdiffx(j,i)=0.0d0 - enddo - enddo do i=1,nfrag_back ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset) c diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index 5707b4b..f27a831 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -92,7 +92,7 @@ C FG slaves receive the WEIGHTS array 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 @@ -307,6 +307,7 @@ C 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 @@ -5345,6 +5346,7 @@ C 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)) @@ -5958,9 +5960,21 @@ c MODELLER restraint function 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' @@ -5971,6 +5985,12 @@ c MODELLER restraint function 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 @@ -5983,16 +6003,26 @@ c MODELLER restraint function 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*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) @@ -6013,18 +6043,22 @@ ccc & distance(i,j,k)**2, "dist(i+1,j+1)=", dist(i+1,j+1), 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)+min_odl - sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist +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 @@ -6033,7 +6067,8 @@ c write(iout,*) "constr_homology=",constr_homology c write(iout,*) i, j, k, "TEST K" enddo - grad_odl3=sum_sgodl/(sum_godl*dij) + grad_odl3=waga_dist*sum_sgodl/(sum_godl*dij) +c grad_odl3=sum_sgodl/(sum_godl*dij) c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2" @@ -6053,41 +6088,67 @@ ccc & ghpbc(jik,i+1), ghpbc(jik,j+1) 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*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 @@ -6097,30 +6158,296 @@ c Gradient 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 + 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 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 + & +waga_d*Erot +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) diff --git a/source/unres/src_MD/gradient_p.F b/source/unres/src_MD/gradient_p.F index 7fec1e8..8606ba4 100644 --- a/source/unres/src_MD/gradient_p.F +++ b/source/unres/src_MD/gradient_p.F @@ -252,10 +252,15 @@ cd enddo C------------------------------------------------------------------------- subroutine cartgrad implicit real*8 (a-h,o-z) +c +c integer iistart,iiend +c include 'DIMENSIONS' + include 'COMMON.LOCAL' #ifdef MPI include 'mpif.h' #endif + include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' @@ -277,30 +282,66 @@ c enddo time00=MPI_Wtime() #endif icg=1 +#ifdef DEBUG +c write (iout,*) "in cartgrad before sum_: duscdiff and duscdiffx" + write (iout,*) "------ Before call sum_ in cargrad ------" + do i=1,nres + write (iout,*) i,(duscdiff(j,i),j=1,3) + write (iout,*) i,(duscdiffx(j,i),j=1,3) +c write (iout,*) "nphi+i",nphi+i," gloc",gloc(nphi+i,icg) +c write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3) +c write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3) + enddo +#endif call sum_gradient +c #ifdef TIMING #endif c do i=1,nres c write (iout,*) "checkgrad", gloc_sc(1,i,icg),gloc(i,icg) c enddo -cd write (iout,*) "After sum_gradient" -cd do i=1,nres-1 -cd write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) -cd write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) -cd enddo +#ifdef DEBUG + write (iout,*) "------ After sum_gradient in cartgrad ------" +c do i=1,nres-1 + do i=1,nres + write (iout,*) "nphi+i",nphi+i," gloc",gloc(nphi+i,icg) + 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 +#endif c If performing constraint dynamics, add the gradients of the constraint energy - if(usampl.and.totT.gt.eq_time) then +#ifdef DEBUG + write (iout,*) "in cartgrad: dutheta, duscdiff and duscdiffx" + do i=1,nres + write (iout,*) i,dutheta(i) + write (iout,*) i,(duscdiff(j,i),j=1,3) + write (iout,*) i,(duscdiffx(j,i),j=1,3) + enddo +#endif +c if(usampl.and.totT.gt.eq_time) then + if(usampl.and.totT.gt.eq_time .or. constr_homology.gt.0) then +c +c Setting suited bounds for HM restrs +c do i=1,nct do j=1,3 gradc(j,i,icg)=gradc(j,i,icg)+dudconst(j,i)+duscdiff(j,i) gradx(j,i,icg)=gradx(j,i,icg)+dudxconst(j,i)+duscdiffx(j,i) enddo +#ifdef DEBUG + write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3) + write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3) +#endif enddo do i=1,nres-3 gloc(i,icg)=gloc(i,icg)+dugamma(i) enddo +c do i=1,nres-2 gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i) +#ifdef DEBUG + write (iout,*) "nphi+i",nphi+i," gloc",gloc(nphi+i,icg) +#endif enddo endif #ifdef TIMING @@ -389,6 +430,17 @@ C enddo enddo enddo +c +c Initialize the gradients of local restraints +c + do i=1,nres + dutheta(i)=0.0d0 + dugamma(i)=0.0d0 + do j=1,3 + duscdiff(j,i)=0.0d0 + duscdiffx(j,i)=0.0d0 + enddo + enddo C C Initialize the gradient of local energy terms. C diff --git a/source/unres/src_MD/initialize_p.F b/source/unres/src_MD/initialize_p.F index d02ebd1..ea446f7 100644 --- a/source/unres/src_MD/initialize_p.F +++ b/source/unres/src_MD/initialize_p.F @@ -37,8 +37,9 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.MINIM' include 'COMMON.DERIV' include 'COMMON.SPLITELE' + include 'COMMON.MD' c Common blocks from the diagonalization routines - COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400) + COMMON /IOFILE/ IR,IW,IPP,IJK,IPK,IDAF,NAV,IODA(400) COMMON /MACHSW/ KDIAG,ICORFL,IXDR logical mask_r c real*8 text1 /'initial_i'/ @@ -245,6 +246,7 @@ C #ifndef SPLITELE nprint_ene=nprint_ene-1 #endif + usampl=.true. return end c------------------------------------------------------------------------- diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index d21d3b9..f1cc85b 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -373,7 +373,6 @@ C 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) @@ -2631,103 +2630,325 @@ c------------------------------------------------------------------------------- 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 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) - - write (iout,*) "nnt",nnt," nct",nct - call flush(iout) + 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 + 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----------------------------------------------------------------------