From: Cezary Czaplewski Date: Fri, 22 May 2020 21:49:39 +0000 (+0200) Subject: Adam's corrections X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=79fd686e9d9aaa23e4fda802ee16379470e961d3 Adam's corrections --- diff --git a/source/cluster/wham/src-HCD-5D/COMMON.CHAIN b/source/cluster/wham/src-HCD-5D/COMMON.CHAIN index 9de64dd..2b481a5 100644 --- a/source/cluster/wham/src-HCD-5D/COMMON.CHAIN +++ b/source/cluster/wham/src-HCD-5D/COMMON.CHAIN @@ -1,6 +1,6 @@ integer nres,nres0,nsup,nstart_sup,nend_sup,nstart_seq, & nchain,chain_border,chain_length,ireschain,npermchain, - & tabpermchain,ishift_pdb,iz_sc + & tabpermchain,ishift_pdb,iz_sc,nres_chomo double precision c,cref,crefjlee,cref_pdb,dc,xloc,xrot,dc_norm, & t,r,prod,rt,chomo common /chain/ c(3,maxres2+2),dc(3,maxres2),xloc(3,maxres), @@ -18,4 +18,5 @@ & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad, & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick - common /chomo_models/ chomo(3,maxres2+2,max_template) + common /chomo_models/ chomo(3,maxres2+2,max_template), + & nres_chomo(max_template) diff --git a/source/cluster/wham/src-HCD-5D/energy_p_new.F b/source/cluster/wham/src-HCD-5D/energy_p_new.F index 27e944b..5d07d5d 100644 --- a/source/cluster/wham/src-HCD-5D/energy_p_new.F +++ b/source/cluster/wham/src-HCD-5D/energy_p_new.F @@ -10203,16 +10203,21 @@ c enddo c min_odl=minval(distancek) - do kk=1,constr_homology - if(l_homo(kk,ii)) then - min_odl=distancek(kk) - exit - endif - enddo - do kk=1,constr_homology - if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + if (nexl.gt.0) then + min_odl=0.0d0 + else + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) & min_odl=distancek(kk) - enddo + enddo + endif + c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij diff --git a/source/cluster/wham/src-HCD-5D/read_constr_homology.F b/source/cluster/wham/src-HCD-5D/read_constr_homology.F index b188deb..6ae3ef4 100644 --- a/source/cluster/wham/src-HCD-5D/read_constr_homology.F +++ b/source/cluster/wham/src-HCD-5D/read_constr_homology.F @@ -1,5 +1,5 @@ subroutine read_constr_homology - + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -27,16 +27,19 @@ 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, ii_in_use(maxdim),i_tmp,idomain_tmp + integer ki,i,ii,ik,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp, + & lim_theta,lim_xx,irec,iistart,iishift,i10,i01 + double precision distal integer idomain(max_template,maxres) - logical lprn /.true./ + logical lfirst integer ilen external ilen logical liiflag + integer nres_temp c c FP - Nov. 2014 Temporary specifications for new vars c - double precision rescore_tmp,x12,y12,z12,rescore2_tmp + double precision rescore_tmp,x12,y12,z12,rescore2_tmp, & rescore3_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 @@ -142,6 +145,7 @@ c tpl_k_rescore="template"//kic2//".sco" unres_pdb=.false. + nres_temp=nres if (read2sigma) then call readpdb_template(k) close(ipdbin) @@ -149,15 +153,16 @@ c call readpdb(out_template_coord) close(ipdbin) endif + nres_chomo(k)=nres + nres=nres_temp -c 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 + do i=1,nres_chomo(k) write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), & (crefjlee(j,i+nres),j=1,3) enddo @@ -170,7 +175,7 @@ c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) - do i=1,2*nres + do i=1,2*nres_chomo(k) do j=1,3 c(j,i)=cref(j,i) c write (iout,*) "c(",j,i,") =",c(j,i) @@ -259,6 +264,8 @@ c & constr_homology enddo lim_odl=ii endif +c write (iout,*) "Distance restraints set" +c call flush(iout) c c Theta, dihedral and SC retraints c @@ -294,11 +301,13 @@ 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 if (sigma_dih(k,i).ne.0) - & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) + & 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 +c write (iout,*) "Dihedral angle restraints set" +c call flush(iout) if (waga_theta.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_theta,status='old') @@ -334,6 +343,8 @@ c rescore(k,i-2) ! right expression ? c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif +c write (iout,*) "Angle restraints set" +c call flush(iout) if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') @@ -359,7 +370,6 @@ 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) -c sigma_d(k,i)=rescore(k,i) ! right expression ? sigma_d(k,i)=rescore3(k,i) ! right expression ? if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) @@ -370,57 +380,70 @@ c read (ientin,*) sigma_d(k,i) ! 1st variant enddo endif enddo +c write (iout,*) "SC restraints set" +c call flush(iout) c c remove distance restraints not used in any model from the list c shift data in all arrays c +c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. + lfirst=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 - if (ii_in_use(ii).eq.0.and.liiflag) then +c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 +c & .and. distal.le.dist2_cut ) then +c write (iout,*) "i",i," j",j," ii",ii +c call flush(iout) + if (ii_in_use(ii).eq.0.and.liiflag.or. + & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then liiflag=.false. - iistart=ii + i10=ii + if (lfirst) then + lfirst=.false. + iistart=ii + else + if(i10.eq.lim_odl) i10=i10+1 + do ki=0,i10-i01-1 + ires_homo(iistart+ki)=ires_homo(ki+i01) + jres_homo(iistart+ki)=jres_homo(ki+i01) + ii_in_use(iistart+ki)=ii_in_use(ki+i01) + do k=1,constr_homology + odl(k,iistart+ki)=odl(k,ki+i01) + sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01) + l_homo(k,iistart+ki)=l_homo(k,ki+i01) + enddo + enddo + iistart=iistart+i10-i01 + endif endif - if (ii_in_use(ii).ne.0.and..not.liiflag.or. - & .not.liiflag.and.ii.eq.lim_odl) then - if (ii.eq.lim_odl) then - iishift=ii-iistart+1 - else - iishift=ii-iistart - endif + if (ii_in_use(ii).ne.0.and..not.liiflag) then + i01=ii liiflag=.true. - do ki=iistart,lim_odl-iishift - ires_homo(ki)=ires_homo(ki+iishift) - jres_homo(ki)=jres_homo(ki+iishift) - ii_in_use(ki)=ii_in_use(ki+iishift) - do k=1,constr_homology - odl(k,ki)=odl(k,ki+iishift) - sigma_odl(k,ki)=sigma_odl(k,ki+iishift) - l_homo(k,ki)=l_homo(k,ki+iishift) - enddo - enddo - ii=ii-iishift - lim_odl=lim_odl-iishift endif enddo enddo + lim_odl=iistart-1 endif - - endif ! .not. klapaucjusz +c write (iout,*) "Removing distances completed" +c call flush(iout) + endif ! .not. klapaucjusz if (constr_homology.gt.0) call homology_partition +c write (iout,*) "After homology_partition" +c call flush(iout) 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 write (iout,*) "After init_int_table" +c call flush(iout) +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 + if (.not.out_template_restr) return cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,*) "Distance restraints from templates" @@ -477,7 +500,9 @@ c---------------------------------------------------------------------- character*500 controlcard integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp integer idomain(max_template,maxres) + integer nres_temp logical lprn /.true./ + logical lfirst integer ilen external ilen logical liiflag @@ -513,7 +538,10 @@ c Read pdb files stop 34 continue unres_pdb=.false. + nres_temp=nres call readpdb_template(k) + nres_chomo(k)=nres + nres=nres_temp c do i=1,2*nres c do j=1,3 c chomo(j,i,k)=c(j,i) @@ -552,6 +580,8 @@ c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) + nres_temp=nres + nres=nres_chomo(k) do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,k) @@ -564,6 +594,7 @@ c write (iout,*) "c(",j,i,") =",c(j,i) thetaref(i)=theta(i) phiref(i)=phi(i) enddo + nres=nres_temp if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 diff --git a/source/cluster/wham/src-HCD-5D/readpdb.F b/source/cluster/wham/src-HCD-5D/readpdb.F index a7bb52b..58c63e4 100644 --- a/source/cluster/wham/src-HCD-5D/readpdb.F +++ b/source/cluster/wham/src-HCD-5D/readpdb.F @@ -98,7 +98,7 @@ c write (2,'(a)') card ! write (iout,*) "Calculating sidechain center iii",iii if (unres_pdb) then do j=1,3 - dc(j,ires+nres)=sccor(j,iii) + dc(j,ires_old)=sccor(j,iii) enddo else call sccenter(ires_old,iii,sccor) @@ -347,17 +347,20 @@ c character*5 atom & ' Gamma' endif endif - do i=1,nres-1 + do i=2,nres iti=itype(i) - if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then +c write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1) + if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and. + & (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then write (iout,'(a,i4)') 'Bad Cartesians for residue',i -ctest stop +c stop endif - vbld(i+1)=dist(i,i+1) - vbld_inv(i+1)=1.0d0/vbld(i+1) - if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1) + vbld(i)=dist(i-1,i) + vbld_inv(i)=1.0d0/vbld(i) + theta(i+1)=alpha(i-1,i,i+1) if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) enddo + c if (unres_pdb) then c if (itype(1).eq.ntyp1) then c theta(3)=90.0d0*deg2rad @@ -806,7 +809,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue endif endif C Calculate internal coordinates. - if (lprn) then + if (out_template_coord) then write (iout,'(/a)') & "Cartesian coordinates of the reference structure" write (iout,'(a,3(3x,a5),5x,3(3x,a5))') @@ -818,6 +821,7 @@ C Calculate internal coordinates. enddo endif C Calculate internal coordinates. +#ifdef DEBUG write (iout,'(a)') & "Backbone and SC coordinates as read from the PDB" do ires=1,nres @@ -825,7 +829,8 @@ C Calculate internal coordinates. & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3), & (c(j,nres+ires),j=1,3) enddo - call int_from_cart(.true.,.false.) +#endif + call int_from_cart(.true.,out_template_coord) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) diff --git a/source/unres/src-HCD-5D/MREMD.F b/source/unres/src-HCD-5D/MREMD.F index 38db8a8..9191402 100644 --- a/source/unres/src-HCD-5D/MREMD.F +++ b/source/unres/src-HCD-5D/MREMD.F @@ -62,6 +62,7 @@ cold integer nup(0:maxprocs),ndown(0:maxprocs) & ene_i_iex,xxx,tmp,econstr_temp_iex,econstr_temp_i integer iran_num double precision ran_number + integer i_econstr/20/ cdeb imin_itime_old=0 ntwx_cache=0 @@ -186,8 +187,9 @@ cd write (*,*) me," After broadcast: file_exist",file_exist do i=nnt,nct if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i))) & *dsqrt(gamsc(iabs(itype(i)))) - enddo + enddo endif + endif if(me.eq.king) then @@ -250,20 +252,20 @@ cd write (*,*) me," After broadcast: file_exist",file_exist write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":", & (ndowna(i,il),i=1,ndowna(0,il)) enddo + stdfp=dsqrt(2*Rb*t_bath/d_time) do i=1,ntyp stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) enddo if (lang.gt.0 .and. .not.surfarea) then - do i=nnt,nct-1 - stdforcp(i)=stdfp*dsqrt(gamp) - enddo + do i=nnt,nct-1 + stdforcp(i)=stdfp*dsqrt(gamp) + enddo do i=nnt,nct if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i))) & *dsqrt(gamsc(iabs(itype(i)))) - enddo - endif - + enddo + endif ELSE IF (.not.(rest.and.file_exist)) THEN do il=1,remd_m(1) ifirst(il)=il @@ -319,7 +321,7 @@ cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep) iset=i2set(me) if(me.eq.king.or..not.out1file) & write(iout,*) me,"iset=",iset,"t_bath=",t_bath - call flush(iout) +c call flush(iout) endif c stdfp=dsqrt(2*Rb*t_bath/d_time) @@ -721,7 +723,7 @@ ctime call flush(iout) write(iout,*) "MIN ii_write=",ii_write endif endif - call flush(iout) +c call flush(iout) endif if (synflag) then c Update the time safety limiy @@ -855,7 +857,7 @@ cd end time05=MPI_WTIME() if (me.eq.king .or. .not. out1file) then write(iout,*) 'REMD writing traj time=',time05-time04 - call flush(iout) +c call flush(iout) endif @@ -1089,12 +1091,12 @@ cd write(iout,*) "i_dir=",i_dir iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1) c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned c on failed replica exchange attempt - econstr_temp_i=remd_ene(20,i) - econstr_temp_iex=remd_ene(20,iex) + econstr_temp_i=remd_ene(i_econstr,i) + econstr_temp_iex=remd_ene(i_econstr,iex) c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials) if (adaptive) then - remd_ene(20,i)=remd_ene(n_ene+5,i) - remd_ene(20,iex)=remd_ene(n_ene+6,iex) + remd_ene(i_econstr,i)=remd_ene(n_ene+5,i) + remd_ene(i_econstr,iex)=remd_ene(n_ene+6,iex) endif elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then @@ -1103,10 +1105,10 @@ c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials) i_iset1=i_iset+1 i_mset1=iran_num(1,mset(i_iset1)) iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1) - econstr_temp_i=remd_ene(20,i) - econstr_temp_iex=remd_ene(20,iex) - remd_ene(20,i)=remd_ene(n_ene+3,i) - remd_ene(20,iex)=remd_ene(n_ene+4,iex) + econstr_temp_i=remd_ene(i_econstr,i) + econstr_temp_iex=remd_ene(i_econstr,iex) + remd_ene(i_econstr,i)=remd_ene(n_ene+3,i) + remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex) elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then @@ -1115,14 +1117,14 @@ c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials) i_iset1=i_iset+1 i_mset1=iran_num(1,mset(i_iset1)) iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1) - econstr_temp_i=remd_ene(20,i) - econstr_temp_iex=remd_ene(20,iex) + econstr_temp_i=remd_ene(i_econstr,i) + econstr_temp_iex=remd_ene(i_econstr,iex) if (adaptive) then - remd_ene(20,i)=remd_ene(n_ene+7,i) - remd_ene(20,iex)=remd_ene(n_ene+8,iex) + remd_ene(i_econstr,i)=remd_ene(n_ene+7,i) + remd_ene(i_econstr,iex)=remd_ene(n_ene+8,iex) else - remd_ene(20,i)=remd_ene(n_ene+3,i) - remd_ene(20,iex)=remd_ene(n_ene+4,iex) + remd_ene(i_econstr,i)=remd_ene(n_ene+3,i) + remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex) endif else goto 444 @@ -1222,8 +1224,8 @@ cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx else remd_ene(0,iex)=ene_iex_iex remd_ene(0,i)=ene_i_i - remd_ene(20,iex)=econstr_temp_iex - remd_ene(20,i)=econstr_temp_i + remd_ene(i_econstr,iex)=econstr_temp_iex + remd_ene(i_econstr,i)=econstr_temp_i endif cd do il=1,nset @@ -1259,7 +1261,7 @@ c------------------------------------- enddo endif - call flush(iout) +c call flush(iout) cd write (iout,'(a6,100i4)') "ifirst", cd & (ifirst(i),i=1,remd_m(1)) @@ -1305,7 +1307,6 @@ co & " rescaling weights with temperature",t_bath do i=1,ntyp stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) enddo - c Compute the standard deviations of stochastic forces for Langevin dynamics c if the friction coefficients do not depend on surface area if (lang.gt.0 .and. .not.surfarea) then @@ -1317,11 +1318,11 @@ c if the friction coefficients do not depend on surface area & *dsqrt(gamsc(iabs(itype(i)))) enddo endif -cde write(iout,*) 'REMD after',me,t_bath +cde write(iout,*) 'REMD after',me,t_bath time08=MPI_WTIME() if (me.eq.king .or. .not. out1file) then write(iout,*) 'REMD exchange time=',time08-time00 - call flush(iout) +c call flush(iout) endif endif enddo diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index 5dd78ea..16c0f37 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -3559,13 +3559,12 @@ c---------------------------------------------------------------------- character*24 model_ki_dist, model_ki_angle character*500 controlcard integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp, - & ik,ll,ii,kk,iistart,lim_xx + & ik,ll,ii,kk,iistart,iishift,lim_xx double precision distal logical lprn /.true./ integer ilen external ilen - logical liiflag,lfirst - integer i01,i10 + logical liiflag c c double precision rescore_tmp,x12,y12,z12,rescore2_tmp @@ -3756,39 +3755,36 @@ c if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. - lfirst=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 - if (ii_in_use(ii).eq.0.and.liiflag.or. - & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then + if (ii_in_use(ii).eq.0.and.liiflag) then liiflag=.false. - i10=ii - if (lfirst) then - lfirst=.false. - iistart=ii - else - if(i10.eq.lim_odl) i10=i10+1 - do ki=0,i10-i01-1 - ires_homo(iistart+ki)=ires_homo(ki+i01) - jres_homo(iistart+ki)=jres_homo(ki+i01) - ii_in_use(iistart+ki)=ii_in_use(ki+i01) - do k=1,constr_homology - odl(k,iistart+ki)=odl(k,ki+i01) - sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01) - l_homo(k,iistart+ki)=l_homo(k,ki+i01) - enddo - enddo - iistart=iistart+i10-i01 - endif + iistart=ii endif - if (ii_in_use(ii).ne.0.and..not.liiflag) then - i01=ii + if (ii_in_use(ii).ne.0.and..not.liiflag.or. + & .not.liiflag.and.ii.eq.lim_odl) then + if (ii.eq.lim_odl) then + iishift=ii-iistart+1 + else + iishift=ii-iistart + endif liiflag=.true. + do ki=iistart,lim_odl-iishift + ires_homo(ki)=ires_homo(ki+iishift) + jres_homo(ki)=jres_homo(ki+iishift) + ii_in_use(ki)=ii_in_use(ki+iishift) + do k=1,constr_homology + odl(k,ki)=odl(k,ki+iishift) + sigma_odl(k,ki)=sigma_odl(k,ki+iishift) + l_homo(k,ki)=l_homo(k,ki+iishift) + enddo + enddo + ii=ii-iishift + lim_odl=lim_odl-iishift endif enddo enddo - lim_odl=iistart-1 endif return diff --git a/source/wham/src-HCD-5D/COMMON.CHAIN b/source/wham/src-HCD-5D/COMMON.CHAIN index dfffc78..7b79a58 100644 --- a/source/wham/src-HCD-5D/COMMON.CHAIN +++ b/source/wham/src-HCD-5D/COMMON.CHAIN @@ -1,6 +1,6 @@ integer nres,nres0,nsup,nstart_sup,nend_sup,nstart_seq, & ishift_pdb,chain_length,chain_border,chain_border1,ichanres, - & tabpermchain,nchain ,npermchain,ireschain,iz_sc + & tabpermchain,nchain ,npermchain,ireschain,iz_sc,nres_chomo double precision c,cref,crefjlee,dc,xloc,xrot,dc_norm,t,r,prod,rt, & rmssing,anatemp,chomo common /chain/ c(3,maxres2+2),dc(3,maxres2),xloc(3,maxres), @@ -18,4 +18,5 @@ & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad, & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick - common /chomo_models/ chomo(3,maxres2+2,max_template) + common /chomo_models/ chomo(3,maxres2+2,max_template), + & nres_chomo(max_template) diff --git a/source/wham/src-HCD-5D/energy_p_new.F b/source/wham/src-HCD-5D/energy_p_new.F index bd69774..6105156 100644 --- a/source/wham/src-HCD-5D/energy_p_new.F +++ b/source/wham/src-HCD-5D/energy_p_new.F @@ -4609,16 +4609,20 @@ c enddo c min_odl=minval(distancek) - do kk=1,constr_homology - if(l_homo(kk,ii)) then - min_odl=distancek(kk) - exit - endif - enddo - do kk=1,constr_homology - if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + if (nexl.gt.0) then + min_odl=0.0d0 + else + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) & min_odl=distancek(kk) - enddo + enddo + endif c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij diff --git a/source/wham/src-HCD-5D/read_constr_homology.F b/source/wham/src-HCD-5D/read_constr_homology.F index dece50f..fa81b80 100644 --- a/source/wham/src-HCD-5D/read_constr_homology.F +++ b/source/wham/src-HCD-5D/read_constr_homology.F @@ -1,5 +1,6 @@ subroutine read_constr_homology + implicit none include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.FREE' @@ -29,12 +30,15 @@ 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, ii_in_use(maxdim),i_tmp,idomain_tmp + integer ki,i,ii,ik,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp, + & lim_theta,lim_xx,irec,iistart,iishift,i10,i01 + double precision distal integer idomain(max_template,maxres) - logical lprn /.true./ + logical lfirst integer ilen external ilen logical liiflag + integer nres_temp c c FP - Nov. 2014 Temporary specifications for new vars c @@ -143,6 +147,7 @@ c tpl_k_rescore="template"//kic2//".sco" unres_pdb=.false. + nres_temp=nres if (read2sigma) then call readpdb_template(k) close(ipdbin) @@ -150,27 +155,29 @@ c call readpdb close(ipdbin) endif + nres_chomo(k)=nres + nres=nres_temp -c call readpdb - do i=1,2*nres + do i=1,2*nres_chomo(k) do j=1,3 crefjlee(j,i)=c(j,i) enddo enddo #ifdef DEBUG - do i=1,nres + do i=1,nres_chmo(k) write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), & (crefjlee(j,i+nres),j=1,3) enddo write (iout,*) "read_constr_homology: after reading pdb file" call flush(iout) #endif + c c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) - do i=1,2*nres + do i=1,2*nres_chomo(k) do j=1,3 c(j,i)=cref(j,i) c write (iout,*) "c(",j,i,") =",c(j,i) @@ -259,6 +266,8 @@ c & constr_homology enddo lim_odl=ii endif +c write (iout,*) "Distance restraints set" +c call flush(iout) c c Theta, dihedral and SC retraints c @@ -294,11 +303,13 @@ 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 if (sigma_dih(k,i).ne.0) - & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) + & 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 +c write (iout,*) "Dihedral angle restraints set" +c call flush(iout) if (waga_theta.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_theta,status='old') @@ -334,6 +345,8 @@ c rescore(k,i-2) ! right expression ? c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif +c write (iout,*) "Angle restraints set" +c call flush(iout) if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') @@ -359,7 +372,6 @@ 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) -c sigma_d(k,i)=rescore(k,i) ! right expression ? sigma_d(k,i)=rescore3(k,i) ! right expression ? if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) @@ -370,57 +382,70 @@ c read (ientin,*) sigma_d(k,i) ! 1st variant enddo endif enddo +c write (iout,*) "SC restraints set" +c call flush(iout) c c remove distance restraints not used in any model from the list c shift data in all arrays c +c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. + lfirst=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 - if (ii_in_use(ii).eq.0.and.liiflag) then +c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 +c & .and. distal.le.dist2_cut ) then +c write (iout,*) "i",i," j",j," ii",ii +c call flush(iout) + if (ii_in_use(ii).eq.0.and.liiflag.or. + & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then liiflag=.false. - iistart=ii + i10=ii + if (lfirst) then + lfirst=.false. + iistart=ii + else + if(i10.eq.lim_odl) i10=i10+1 + do ki=0,i10-i01-1 + ires_homo(iistart+ki)=ires_homo(ki+i01) + jres_homo(iistart+ki)=jres_homo(ki+i01) + ii_in_use(iistart+ki)=ii_in_use(ki+i01) + do k=1,constr_homology + odl(k,iistart+ki)=odl(k,ki+i01) + sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01) + l_homo(k,iistart+ki)=l_homo(k,ki+i01) + enddo + enddo + iistart=iistart+i10-i01 + endif endif - if (ii_in_use(ii).ne.0.and..not.liiflag.or. - & .not.liiflag.and.ii.eq.lim_odl) then - if (ii.eq.lim_odl) then - iishift=ii-iistart+1 - else - iishift=ii-iistart - endif + if (ii_in_use(ii).ne.0.and..not.liiflag) then + i01=ii liiflag=.true. - do ki=iistart,lim_odl-iishift - ires_homo(ki)=ires_homo(ki+iishift) - jres_homo(ki)=jres_homo(ki+iishift) - ii_in_use(ki)=ii_in_use(ki+iishift) - do k=1,constr_homology - odl(k,ki)=odl(k,ki+iishift) - sigma_odl(k,ki)=sigma_odl(k,ki+iishift) - l_homo(k,ki)=l_homo(k,ki+iishift) - enddo - enddo - ii=ii-iishift - lim_odl=lim_odl-iishift endif enddo enddo + lim_odl=iistart-1 endif - - endif ! .not. klapaucjusz +c write (iout,*) "Removing distances completed" +c call flush(iout) + endif ! .not. klapaucjusz if (constr_homology.gt.0) call homology_partition +c write (iout,*) "After homology_partition" +c call flush(iout) 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 write (iout,*) "After init_int_table" +c call flush(iout) +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 + if (.not.out_template_restr) return cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,*) "Distance restraints from templates" @@ -454,7 +479,7 @@ c ----------------------------------------------------------------- end c---------------------------------------------------------------------- subroutine read_klapaucjusz - + implicit none include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.FREE' @@ -477,8 +502,11 @@ c---------------------------------------------------------------------- character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard - integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp + integer i,ii,ik,ki,j,k,l,ll,ii_in_use(maxdim),i_tmp,idomain_tmp, + & nclust,iistart,iishift,lim_xx + integer nres_temp integer idomain(max_template,maxres) + double precision distal logical lprn /.true./ integer ilen external ilen @@ -498,6 +526,7 @@ c call getenv("FRAGFILE",fragfile) open(ientin,file=fragfile,status="old",err=10) read(ientin,*) constr_homology,nclust + lim_xx=0 l_homo = .false. sigma_theta=0.0 sigma_d=0.0 @@ -514,7 +543,10 @@ c Read pdb files stop 34 continue unres_pdb=.false. + nres_temp=nres call readpdb_template(k) + nres_chomo(k)=nres + nres=nres_temp c do i=1,2*nres c do j=1,3 c chomo(j,i,k)=c(j,i) @@ -535,6 +567,7 @@ c c Loop over clusters c do l=1,nclust +c write (iout,*) "CLUSTER",l," members:",ninclust(l) do ll = 1,ninclust(l) k = inclust(ll,l) @@ -553,18 +586,27 @@ c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) +c write (iout,*) "k",k + nres_temp=nres + nres=nres_chomo(k) do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,k) -c write (iout,*) "c(",j,i,") =",c(j,i) +c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo +c write (iout,*) "read_klapaucjusz: calling int_from_cart" +c call cartprint +c write (iout,*) "idomain" +c write (iout,'(2i5)') (i,idomain(k,i),i=1,nres) call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) +c write (iout,*) "en" do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo + nres=nres_temp if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 diff --git a/source/wham/src-HCD-5D/readpdb.F b/source/wham/src-HCD-5D/readpdb.F index 83c63ab..6f4ba5f 100644 --- a/source/wham/src-HCD-5D/readpdb.F +++ b/source/wham/src-HCD-5D/readpdb.F @@ -27,10 +27,9 @@ C geometry. ibeg=1 ishift1=0 sccalc=.false. - ires=0 do read (ipdbin,'(a80)',end=10) card -c write (iout,'(a)') card +! write (iout,'(a)') card if (card(:5).eq.'HELIX') then nhfrag=nhfrag+1 lsecondary=.true. @@ -59,7 +58,7 @@ c write (iout,'(a)') card iterter(ires_old)=1 ishift1=ishift1+1 ibeg=2 - write (iout,*) "Chain ended",ires,ishift,ires_old +! write (iout,*) "Chain ended",ires,ishift,ires_old if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -92,7 +91,7 @@ c write (2,'(a)') card ! write (iout,*) "Calculating sidechain center iii",iii if (unres_pdb) then do j=1,3 - dc(j,ires+nres)=sccor(j,iii) + dc(j,ires_old)=sccor(j,iii) enddo else call sccenter(ires_old,iii,sccor) @@ -120,7 +119,7 @@ c write (iout,*) "BEG ires",ires ! Start a new chain ishift=-ires_old+ires-1 !!!!! ishift1=ishift1-1 !!!!! - write (iout,*) "New chain started",ires,ishift,ishift1,"!" +! write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 ires_old=ires ibeg=0 @@ -326,7 +325,7 @@ c write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1) if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and. & (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then write (iout,'(a,i4)') 'Bad Cartesians for residue',i - stop +c stop endif vbld(i)=dist(i-1,i) vbld_inv(i)=1.0d0/vbld(i) @@ -805,8 +804,8 @@ C Calculate internal coordinates. endif C Calculate internal coordinates. c call int_from_cart1(.false.) - call int_from_cart(.true.,.true.) - call sc_loc_geom(.true.) + call int_from_cart(.true.,out_template_coord) + call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i)