X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=73ffe83c721c92f3590f0983cf572d4e3959f5b2;hb=fb112c2c0d337f5433d28f2c66ecf420102e63f0;hp=671943243338236939dae878b8c9e372b89406a1;hpb=b3571ff407666d8bc4f750cbd8384322cd1d7b61;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 6719432..73ffe83 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -8,6 +8,7 @@ include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' logical file_exist C Read force-field parameters except weights call parmread @@ -130,6 +131,7 @@ C Set up the time limit (caution! The time must be input in minutes!) 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) @@ -373,7 +375,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) @@ -387,6 +388,11 @@ C 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 @@ -427,6 +433,10 @@ c if performing umbrella sampling, fragments constrained are read from the frag 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 @@ -628,6 +638,7 @@ C common /pizda/ itype_pdb logical seq_comp,fail double precision energia(0:n_ene) + integer ilen external ilen C @@ -923,6 +934,17 @@ C 12/1/95 Added weight for the multi-body term WCORR 34 continue c print *,'Begin reading pdb data' call readpdb + 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 c print *,'Finished reading pdb data' if(me.eq.king.or..not.out1file) & write (iout,'(a,i3,a,i3)')'nsup=',nsup, @@ -1065,6 +1087,7 @@ C Juyong:READ read_info 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 @@ -1072,6 +1095,7 @@ C both routines should be in dfa.F file!! call read_dfa_info print*, 'read_dfa_info finished!' endif +#endif C if (pdbref) then if(me.eq.king.or..not.out1file) @@ -1166,6 +1190,44 @@ c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup if (constr_homology.gt.0) then 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)=crefjlee(j,i) + 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 @@ -1183,6 +1245,8 @@ C initial geometry. 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) @@ -1201,7 +1265,7 @@ C initial geometry. enddo endif enddo - return +c return else call read_angles(inp,*36) endif @@ -1285,6 +1349,8 @@ C Generate distance constraints, if the PDB structure is to be regularized. 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 @@ -2263,6 +2329,8 @@ c print *,"Processor",myrank," fg_rank",fg_rank & //'.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) @@ -2281,6 +2349,7 @@ c print *,"Processor",myrank," fg_rank",fg_rank 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) @@ -2441,33 +2510,36 @@ c------------------------------------------------------------------------------- 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 @@ -2634,118 +2706,413 @@ 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 - character*3200 controlcard1 - integer ki, i, j, k, l + integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp logical lprn /.true./ - +c +c FP - Nov. 2014 Temporary specifications for new vars +c + double precision rescore_tmp,x12,y12,z12,rescore2_tmp + double precision, dimension (max_template,maxres) :: rescore + double precision, dimension (max_template,maxres) :: rescore2 + 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) + read2sigma=(index(controlcard,'READ2SIGMA').gt.0) if (homol_nset.gt.1)then call card_concat(controlcard) - read(controlcard,*) (waga_dist(i),i=1,homol_nset) - call card_concat(controlcard) - read(controlcard,*) (waga_angle(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) - enddo + 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 - call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0) - call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0) + waga_homology(1)=1.0 endif - write (iout,*) "nnt",nnt," nct",nct - call flush(iout) +cd write (iout,*) "nnt",nnt," nct",nct +cd 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 - enddo +c +c New +c + lim_theta=0 + lim_xx=0 +c + write(iout,*) 'nnt=',nnt,'nct=',nct +c + do i = nnt,nct + do k=1,constr_homology + idomain(k,i)=0 enddo enddo - do i=1,nres-3 - do ki=1,constr_homology - dih(ki,i)=0.0d0 - sigma_dih(ki,i)=0.0d0 + + ii=0 + do i = nnt,nct-2 + do j=i+2,nct + ii=ii+1 + ii_in_use(ii)=0 enddo enddo - do ki=1,constr_homology - write(kic2,'(i2)') ki - if (ki.le.9) kic2="0"//kic2(2: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 + + do k=1,constr_homology + +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 +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') + if (nnt.gt.1) rescore(k,1)=0.0d0 + do irec=nnt,maxdim ! loop for reading res sim + if (read2sigma) then + read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, + & idomain_tmp + i_tmp=i_tmp+nnt-1 + idomain(k,i_tmp)=idomain_tmp + rescore(k,i_tmp)=0.5d0*(rescore_tmp+0.5d0) + rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0) + else + idomain(k,irec)=1 + 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) + endif + enddo + 1401 continue + close (ientin) + + 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' + + unres_pdb=.false. + if (read2sigma) then + call readpdb_template(k) + else + call readpdb + endif +c +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 - 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 + + + + 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 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 + + 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 + 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 (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) ! 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 + endif + 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 sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)* +c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2) + +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 + else + ii=ii+1 + l_homo(k,ii)=.false. + endif + enddo + enddo + lim_odl=ii + 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. + if (idomain(k,i).eq.0) then + sigma_dih(k,i)=0.0 + cycle + endif + 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))/4.0 +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) + enddo + lim_dih=nct-nnt-2 + 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. + if (idomain(k,i).eq.0) then + sigma_theta(k,i)=0.0 + cycle + endif + 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)/3.0 + 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) enddo -c write (iout,*) "ki",ki - if (ki.gt.constr_homology) cycle + endif + lim_theta=nct-nnt-1 + + 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 + + + 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) 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) +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 + enddo + lim_xx=nct-nnt+1 + endif + enddo +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 - 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 + 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 - enddo - lim_odl=ii + 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, +cd & "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,*) "Distance restraints from templates" - do ii=1,lim_odl - write(iout,'(3i5,10(2f8.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 +cd 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,'(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 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) - - + enddo + 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 + endif +c ----------------------------------------------------------------- return end c----------------------------------------------------------------------