From 3d6f9e7811a3341955e6fdda721ab6de34efe252 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Fri, 4 Mar 2016 17:57:37 +0900 Subject: [PATCH] Adam's changes to wham and cluster following previous commit XG.. and ..GX sequence is automatically corrected also in cluster_wham it is possible to select number of clusters instead of cutoff --- source/cluster/wham/src/COMMON.CLUSTER | 4 +-- source/cluster/wham/src/COMMON.CONTROL | 9 ++++-- source/cluster/wham/src/energy_p_new.F | 2 -- source/cluster/wham/src/main_clust.F | 48 ++++++++++++++++++++++++-------- source/cluster/wham/src/probabl.F | 11 +++++--- source/cluster/wham/src/readpdb.F | 2 +- source/cluster/wham/src/readrtns.F | 36 ++++++++++++++++++------ source/cluster/wham/src/wrtclust.f | 11 ++++++-- source/wham/src/energy_p_new.F | 2 -- source/wham/src/molread_zs.F | 13 +++++++++ 10 files changed, 103 insertions(+), 35 deletions(-) diff --git a/source/cluster/wham/src/COMMON.CLUSTER b/source/cluster/wham/src/COMMON.CLUSTER index 4477d19..f1ad0fd 100644 --- a/source/cluster/wham/src/COMMON.CLUSTER +++ b/source/cluster/wham/src/COMMON.CLUSTER @@ -4,11 +4,11 @@ real*4 diss,allcart double precision enetb,entfac,totfree,energy,rmstb integer ncut,ngr,licz,nconf,iass,icc,mult,list_conf, - & nss_all,ihpb_all,jhpb_all,iass_tot,iscore,nprop + & nss_all,ihpb_all,jhpb_all,iass_tot,iscore,nprop,nclust common /clu/ diss(maxdist),energy(0:maxconf), & enetb(max_ene,maxstr_proc),ecut, & entfac(maxconf),totfree(0:maxconf),totfree_gr(maxgr), - & rcutoff(max_cut+1),ncut,min_var,tree,plot_tree,lgrp + & rcutoff(max_cut+1),ncut,nclust,min_var,tree,plot_tree,lgrp common /clu1/ ngr,licz(maxgr),nconf(maxgr,maxingr),iass(maxgr), & iass_tot(maxgr,max_cut),list_conf(maxconf) common /alles/ allcart(3,maxres2,maxstr_proc),rmstb(maxconf), diff --git a/source/cluster/wham/src/COMMON.CONTROL b/source/cluster/wham/src/COMMON.CONTROL index 2b9b9e5..fe7c63a 100644 --- a/source/cluster/wham/src/COMMON.CONTROL +++ b/source/cluster/wham/src/COMMON.CONTROL @@ -5,11 +5,14 @@ real*8 waga_dist, waga_angle,waga_theta, waga_d, dist_cut logical refstr,pdbref,punch_dist,print_dist,caonly,lside, & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx, - & with_dihed_constr,out1file + & with_dihed_constr,out1file,print_homology_restraints, + & print_contact_map,print_homology_models common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2, & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int, & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,constr_dist, - & with_dihed_constr, constr_homology,homol_nset,out1file + & with_dihed_constr, constr_homology,homol_nset,out1file, + & print_contact_map common /cntrlr/ waga_homology(10), - & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,iset + & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,iset, + & print_homology_restraints,print_homology_models diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F index a5ba7fb..8ce7e5b 100644 --- a/source/cluster/wham/src/energy_p_new.F +++ b/source/cluster/wham/src/energy_p_new.F @@ -5166,11 +5166,9 @@ c 3 = SC...Ca...Ca...SCi cosphi=dcos(j*tauangle(intertyp,i)) sinphi=dsin(j*tauangle(intertyp,i)) esccor=esccor+v1ij*cosphi+v2ij*sinphi -#define DEBUG #ifdef DEBUG esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi #endif -#undef DEBUG gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci diff --git a/source/cluster/wham/src/main_clust.F b/source/cluster/wham/src/main_clust.F index c31847f..980871c 100644 --- a/source/cluster/wham/src/main_clust.F +++ b/source/cluster/wham/src/main_clust.F @@ -34,7 +34,7 @@ C DIMENSION NN(maxconf),DISNN(maxconf) LOGICAL FLAG(maxconf) integer i,j,k,l,m,n,len,lev,idum,ii,ind,ioffset,jj,icut,ncon, - & it,ncon_work,ind1,ilen + & it,ncon_work,ind1,ilen,is,ie double precision t1,t2,tcpu,difconf real diss_(maxdist) @@ -76,6 +76,12 @@ c if (refstr) call read_ref_structure(*30) print *,'MAIN: nnt=',nnt,' nct=',nct + if (nclust.gt.0) then + PRINTANG(1)=.TRUE. + PRINTPDB(1)=outpdb + printmol2(1)=outmol2 + ncut=0 + else DO I=1,NCUT PRINTANG(I)=.FALSE. PRINTPDB(I)=0 @@ -87,12 +93,21 @@ c if (refstr) call read_ref_structure(*30) printmol2(i)=outmol2 ENDIF ENDDO + endif + if (ncut.gt.0) then write (iout,*) 'Number of cutoffs:',NCUT write (iout,*) 'Cutoff values:' DO ICUT=1,NCUT WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT), & printpdb(icut),printmol2(icut) ENDDO + else if (nclust.gt.0) then + write (iout,'("Number of clusters requested",i5)') nclust + else + if (me.eq.Master) + & write (iout,*) "ERROR: Either nclust or ncut must be >0" + stop + endif DO I=1,NRES-3 MULT(I)=1 ENDDO @@ -247,29 +262,39 @@ C CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT) c CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL) +c 3/3/16 AL: added explicit number of cluters + if (nclust.gt.0) then + is=nclust-1 + ie=nclust-1 + icut=1 + else + is=1 + ie=lev-1 + endif do i=1,maxgr licz(i)=0 enddo icut=1 - i=1 - NGR=i+1 + i=is + NGR=is+1 do j=1,n licz(iclass(j,i))=licz(iclass(j,i))+1 nconf(iclass(j,i),licz(iclass(j,i)))=j c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)), c & nconf(iclass(j,i),licz(iclass(j,i))) enddo - do i=1,lev-1 - +c do i=1,lev-1 + do i=is,ie idum=lev-i DO L=1,LEV IF (HEIGHT(L).EQ.IDUM) GOTO 190 ENDDO 190 IDUM=L - write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM), - & " icut",icut," cutoff",rcutoff(icut) - IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN - WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut) +c write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM), +c & " icut",icut," cutoff",rcutoff(icut) + IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN + if (nclust.le.0) + & WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut) write (iout,'(a,f8.2)') 'Maximum distance found:', & CRITVAL(IDUM) CALL SRTCLUST(ICUT,ncon_work,iT) @@ -282,9 +307,10 @@ c & nconf(iclass(j,i),licz(iclass(j,i))) do l=1,maxgr licz(l)=0 enddo + ii=i-is+1 do j=1,n - licz(iclass(j,i))=licz(iclass(j,i))+1 - nconf(iclass(j,i),licz(iclass(j,i)))=j + licz(iclass(j,ii))=licz(iclass(j,ii))+1 + nconf(iclass(j,ii),licz(iclass(j,ii)))=j c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)), c & nconf(iclass(j,i),licz(iclass(j,i))) cd print *,j,iclass(j,i), diff --git a/source/cluster/wham/src/probabl.F b/source/cluster/wham/src/probabl.F index 2f34f61..62da77f 100644 --- a/source/cluster/wham/src/probabl.F +++ b/source/cluster/wham/src/probabl.F @@ -197,16 +197,20 @@ cc if (wcorr6.eq.0) ecorr6=0.0d0 #ifdef MPI Fdimless_(i)=beta_h(ib)*etot+entfac(ii) totfree_(i)=etot +#ifdef DEBUG + write (iout,*) "etrop", i,ii,ib, + & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i), + & entfac(ii),Fdimless_(i) +#endif #else Fdimless(i)=beta_h(ib)*etot+entfac(ii) totfree(i)=etot -#endif #ifdef DEBUG - write (iout,*) "etrop", i,ii,ib, & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i), & entfac(ii),Fdimless(i) #endif +#endif enddo ! i #ifdef MPI call MPI_Gatherv(Fdimless_(1),scount(me), @@ -226,6 +230,7 @@ cc if (wcorr6.eq.0) ecorr6=0.0d0 entfac(i)=entfac_(i) enddo #endif + #ifdef DEBUG write (iout,*) "The FDIMLESS array before sorting" do i=1,ncon @@ -251,14 +256,12 @@ c write (iout,*) "qfree",qfree sumprob=0.0 do i=1,min0(ncon,maxstr_proc)-1 sumprob=sumprob+dexp(dble(-fdimless(i)+fdimless(1)))/qfree -c#define DEBUG #ifdef DEBUG write (iout,*) 'i=',i,ib,beta_h(ib), & 1.0d0/(1.987d-3*beta_h(ib)),list_conf(i), & totfree(list_conf(i)), & -entfac(list_conf(i)),fdimless(i),sumprob #endif -c#undef DEBUG if (sumprob.gt.prob_limit) goto 122 c if (sumprob.gt.1.00d0) goto 122 nlist=nlist+1 diff --git a/source/cluster/wham/src/readpdb.F b/source/cluster/wham/src/readpdb.F index 8351702..604b840 100644 --- a/source/cluster/wham/src/readpdb.F +++ b/source/cluster/wham/src/readpdb.F @@ -223,7 +223,7 @@ C Calculate internal coordinates. enddo endif C Calculate internal coordinates. - if(me.eq.king.or..not.out1file)then + if(print_homology_models.and.(me.eq.king.or..not.out1file))then write (iout,'(a)') & "Backbone and SC coordinates as read from the PDB" do ires=1,nres diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F index be830b5..037eb1f 100644 --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@ -37,7 +37,8 @@ C min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) - call readi(controlcard,'NCUT',ncut,1) + call readi(controlcard,'NCUT',ncut,0) + call readi(controlcard,'NCLUST',nclust,5) call readi(controlcard,'NSTART',nstart,0) call readi(controlcard,'NEND',nend,0) call reada(controlcard,'ECUT',ecut,10.0d0) @@ -46,7 +47,8 @@ C lgrp=(index(controlcard,'LGRP').gt.0) caonly=(index(controlcard,'CA_ONLY').gt.0) print_dist=(index(controlcard,'PRINT_DIST').gt.0) - call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) + if (ncut.gt.0) + & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) call readi(controlcard,'IOPT',iopt,2) lside = index(controlcard,"SIDE").gt.0 efree = index(controlcard,"EFREE").gt.0 @@ -69,6 +71,11 @@ C call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) write (iout,*) "with_homology_constr ",with_dihed_constr, & " CONSTR_HOMOLOGY",constr_homology + print_homology_restraints= + & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0 + print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0 + print_homology_models= + & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0 #ifdef AIX call flush_(iout) @@ -268,6 +275,19 @@ C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo + if (itype(2).eq.10) then + write (iout,*) + & "Glycine is the first full residue, initial dummy deleted" + do i=1,nres + itype(i)=itype(i+1) + enddo + nres=nres-1 + endif + if (itype(nres).eq.10) then + write (iout,*) + & "Glycine is the last full residue, terminal dummy deleted" + nres=nres-1 + endif print *,nres print '(20i4)',(itype(i),i=1,nres) @@ -362,7 +382,7 @@ C both routines should be in dfa.F file!! endif if (constr_homology.gt.0) then - call read_constr_homology + call read_constr_homology(print_homology_restraints) endif c if (pdbref) then @@ -432,7 +452,7 @@ c endif enddo enddo endif - call contact(.true.,ncont_ref,icont_ref) + call contact(print_contact_map,ncont_ref,icont_ref) endif c Read distance restraints if (constr_dist.gt.0) then @@ -885,7 +905,7 @@ c call flush(iout) end c====------------------------------------------------------------------- - subroutine read_constr_homology + subroutine read_constr_homology(lprn) include 'DIMENSIONS' #ifdef MPI @@ -913,7 +933,7 @@ c & sigma_odl_temp(maxres,maxres,max_template) character*24 model_ki_dist, model_ki_angle character*500 controlcard integer ki, i, j, k, l - logical lprn /.true./ + logical lprn logical unres_pdb c c FP - Nov. 2014 Temporary specifications for new vars @@ -991,7 +1011,7 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d do k=1,constr_homology read(inp,'(a)') pdbfile - write (iout,*) "k ",k," pdbfile ",pdbfile +c write (iout,*) "k ",k," pdbfile ",pdbfile c Next stament causes error upon compilation (?) c if(me.eq.king.or. .not. out1file) c write (iout,'(2a)') 'PDB data will be read from file ', @@ -1026,7 +1046,6 @@ c tpl_k_sigma_d="template"//kic2//".sigma_d" write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), & (crefjlee(j,i+nres),j=1,3) enddo -#endif write (iout,*) "READ HOMOLOGY INFO" write (iout,*) "read_constr_homology x: after reading pdb file" write (iout,*) "waga_homology(",iset,")",waga_homology(iset) @@ -1035,6 +1054,7 @@ c tpl_k_sigma_d="template"//kic2//".sigma_d" write (iout,*) "waga_theta",waga_theta write (iout,*) "waga_d",waga_d write (iout,*) "dist_cut",dist_cut +#endif #ifdef AIX call flush_(iout) #else diff --git a/source/cluster/wham/src/wrtclust.f b/source/cluster/wham/src/wrtclust.f index 6012b05..ec43c77 100644 --- a/source/cluster/wham/src/wrtclust.f +++ b/source/cluster/wham/src/wrtclust.f @@ -84,12 +84,16 @@ C 12/8/93 Estimation of "diameters" of the subsequent families. ave_dim=0.0 amax_dim=0.0 c write (iout,*) "ecut",ecut + emin=totfree(nconf(igr,1)) +c write (2,*) "emin",emin," ecut",ecut do i=2,licz(igr) ii=nconf(igr,i) +c write (2,*) " igr",igr," i",i," ii",ii," totfree",totfree(ii), +c & " emin",emin," diff",totfree(ii)-emin," ecut",ecut if (totfree(ii)-emin .gt. ecut) goto 10 do j=1,i-1 jj=nconf(igr,j) - if (jj.eq.1) exit +c if (jj.eq.1) exit if (ii.lt.jj) then ind=ioffset(ncon,ii,jj) else @@ -112,9 +116,12 @@ c & list_conf(jj),curr_dist & '; average distance in the family:',ave_dim rmsave(igr)=0.0d0 qpart=0.0d0 + emin=totfree(nconf(igr,1)) do i=1,licz(igr) icon=nconf(igr,i) - boltz=dexp(-totfree(icon)) + boltz=dexp(-totfree(icon)+emin) +c write (2,*) "igr",igr," i",i," icon",icon," totfree", +c & totfree(icon)," emin",emin," boltz",boltz," rms",rmstb(icon) rmsave(igr)=rmsave(igr)+boltz*rmstb(icon) qpart=qpart+boltz enddo diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F index 8a1f9c7..81abea5 100644 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@ -5261,11 +5261,9 @@ c 3 = SC...Ca...Ca...SCi cosphi=dcos(j*tauangle(intertyp,i)) sinphi=dsin(j*tauangle(intertyp,i)) esccor=esccor+v1ij*cosphi+v2ij*sinphi -#define DEBUG #ifdef DEBUG esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi #endif -#undef DEBUG gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci diff --git a/source/wham/src/molread_zs.F b/source/wham/src/molread_zs.F index a98f47c..9c999cc 100644 --- a/source/wham/src/molread_zs.F +++ b/source/wham/src/molread_zs.F @@ -57,6 +57,19 @@ C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo + if (itype(2).eq.10) then + write (iout,*) + & "Glycine is the first full residue, initial dummy deleted" + do i=1,nres + itype(i)=itype(i+1) + enddo + nres=nres-1 + endif + if (itype(nres).eq.10) then + write (iout,*) + & "Glycine is the last full residue, terminal dummy deleted" + nres=nres-1 + endif write (iout,*) "Numeric code:" write (iout,'(20i4)') (itype(i),i=1,nres) do i=1,nres-1 -- 1.7.9.5