From 809cbbbaf17b645d881fe5c61387d09109acdf19 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Tue, 16 Feb 2021 14:46:07 +0100 Subject: [PATCH] Adam's cluster wham update --- source/cluster/wham/src-HCD/COMMON.CHAIN | 2 +- source/cluster/wham/src-HCD/COMMON.DFA | 2 +- source/cluster/wham/src-HCD/COMMON.SHIELD | 9 +- source/cluster/wham/src-HCD/chain_symmetry.F | 3 + source/cluster/wham/src-HCD/energy_p_new.F | 2 + source/cluster/wham/src-HCD/geomout.F | 14 ++- source/cluster/wham/src-HCD/read_constr_homology.F | 4 +- source/cluster/wham/src-HCD/read_coords.F | 98 ++++++++++++-------- source/cluster/wham/src-HCD/readrtns.F | 1 + source/cluster/wham/src-HCD/ssMD.F | 6 +- source/cluster/wham/src-HCD/wrtclust.f | 26 +++++- 11 files changed, 114 insertions(+), 53 deletions(-) diff --git a/source/cluster/wham/src-HCD/COMMON.CHAIN b/source/cluster/wham/src-HCD/COMMON.CHAIN index 2b481a5..a17e113 100644 --- a/source/cluster/wham/src-HCD/COMMON.CHAIN +++ b/source/cluster/wham/src-HCD/COMMON.CHAIN @@ -3,7 +3,7 @@ & 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), + common /chain/ c(3,maxres2+2),dc(3,0:maxres2),xloc(3,maxres), & xrot(3,maxres),dc_norm(3,maxres2),nres,nres0 common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres), & rt(3,3,maxres) diff --git a/source/cluster/wham/src-HCD/COMMON.DFA b/source/cluster/wham/src-HCD/COMMON.DFA index 064a7ce..782e8c4 100644 --- a/source/cluster/wham/src-HCD/COMMON.DFA +++ b/source/cluster/wham/src-HCD/COMMON.DFA @@ -11,7 +11,7 @@ C Total : ~ 11 * Nres restraints C C INTEGER IDFAMAX,IDFAMX2,IDFACMD,IDMAXMIN, MAXN - PARAMETER(IDFAMAX=10000,IDFAMX2=1000,IDFACMD=500,IDMAXMIN=500) + PARAMETER(IDFAMAX=25000,IDFAMX2=1000,IDFACMD=500,IDMAXMIN=500) PARAMETER(MAXN=4) real*8 wwdist,wwangle,wwnei parameter(wwdist=1.0d0,wwangle=1.0d0,wwnei=1.0d0) diff --git a/source/cluster/wham/src-HCD/COMMON.SHIELD b/source/cluster/wham/src-HCD/COMMON.SHIELD index 1f96c94..8d89f0b 100644 --- a/source/cluster/wham/src-HCD/COMMON.SHIELD +++ b/source/cluster/wham/src-HCD/COMMON.SHIELD @@ -5,10 +5,11 @@ common /shield/ VSolvSphere,VSolvSphere_div,buff_shield, & long_r_sidechain(ntyp), & short_r_sidechain(ntyp),fac_shield(maxres),wshield, - & grad_shield_side(3,maxcont,-1:maxres),grad_shield(3,-1:maxres), - & grad_shield_loc(3,maxcont,-1:maxres), - & ishield_list(maxres),shield_list(maxcont,maxres), - & ees0plist(maxcont,maxres) + & grad_shield_side(3,maxint_res,-1:maxres), + & grad_shield(3,-1:maxres), + & grad_shield_loc(3,maxint_res,-1:maxres), + & ishield_list(maxres),shield_list(maxint_res,maxres), + & ees0plist(maxint_res,maxres) diff --git a/source/cluster/wham/src-HCD/chain_symmetry.F b/source/cluster/wham/src-HCD/chain_symmetry.F index 1406d1d..8c36855 100644 --- a/source/cluster/wham/src-HCD/chain_symmetry.F +++ b/source/cluster/wham/src-HCD/chain_symmetry.F @@ -7,6 +7,7 @@ c implicit none include "DIMENSIONS" include "COMMON.IOUNITS" + include "COMMON.CONTROL" integer nchain,nres,itype(nres),chain_border(2,maxchain), & chain_length(nchain),itemp(maxchain), & npermchain,tabpermchain(maxchain,maxperm), @@ -42,6 +43,7 @@ c nchain_group=nchain_group+1 iieq=1 iequiv(iieq,nchain_group)=i + if (symetr.eq.1) then do j=i+1,nchain if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle c k=0 @@ -57,6 +59,7 @@ c k=k+1 iieq=iieq+1 iequiv(iieq,nchain_group)=j enddo + endif nequiv(nchain_group)=iieq enddo write(iout,*) "Number of equivalent chain groups:",nchain_group diff --git a/source/cluster/wham/src-HCD/energy_p_new.F b/source/cluster/wham/src-HCD/energy_p_new.F index 4fa79c5..e969ea3 100644 --- a/source/cluster/wham/src-HCD/energy_p_new.F +++ b/source/cluster/wham/src-HCD/energy_p_new.F @@ -3544,6 +3544,8 @@ C Third- and fourth-order contributions from turns common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij j=i+3 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C diff --git a/source/cluster/wham/src-HCD/geomout.F b/source/cluster/wham/src-HCD/geomout.F index be43686..d513764 100644 --- a/source/cluster/wham/src-HCD/geomout.F +++ b/source/cluster/wham/src-HCD/geomout.F @@ -10,10 +10,22 @@ include 'COMMON.SBRIDGE' include 'COMMON.TEMPFAC' character*50 tytul - character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/ dimension ica(maxres) + character*1 chainid(32) /'A','B','C','D','E','F','G','H','I','J', + & 'K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z', + & '1','2','3','4','5','6'/ write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20), & ' ENERGY ',etot,' RMS ',rmsd + do i=1,nss + ici=ireschain(ihpb(i)-nres) + icj=ireschain(jhpb(i)-nres) + iri=ihpb(i)-chain_border(1,ici)+1-nres + irj=jhpb(i)-chain_border(1,icj)+1-nres +c write (iout,*) ihpb(i),ici,iri,jhpb(i),icj,irj + write(ipdb,'(a6,i4,1x,a3,1x,a1,i5,4x,a3,1x,a1,i5,38x,f5.2)') + & 'SSBOND',i,'CYS',chainid(ici),iri,'CYS',chainid(icj),irj, + & dist(ihpb(i),jhpb(i)) + enddo iatom=0 ichain=1 ires=0 diff --git a/source/cluster/wham/src-HCD/read_constr_homology.F b/source/cluster/wham/src-HCD/read_constr_homology.F index 0b265fa..f2fb1d0 100644 --- a/source/cluster/wham/src-HCD/read_constr_homology.F +++ b/source/cluster/wham/src-HCD/read_constr_homology.F @@ -257,8 +257,8 @@ c & constr_homology endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else - ii=ii+1 - l_homo(k,ii)=.false. +c ii=ii+1 +c l_homo(k,ii)=.false. endif enddo enddo diff --git a/source/cluster/wham/src-HCD/read_coords.F b/source/cluster/wham/src-HCD/read_coords.F index facbc27..20abce5 100644 --- a/source/cluster/wham/src-HCD/read_coords.F +++ b/source/cluster/wham/src-HCD/read_coords.F @@ -49,8 +49,8 @@ c Set the scratchfile names #endif c 1/27/05 AL Change stored coordinates to single precision and don't store c energy components in the binary databases. - lenrec=12*(nres+nct-nnt+1)+4*(2*nss+2)+16 - lenrec_in=12*(nres+nct-nnt+1)+4*(2*nss+2)+24 + lenrec=12*(nres+nct-nnt+1)+4*(ns+2)+16 + lenrec_in=12*(nres+nct-nnt+1)+4*(ns+2)+24 #ifdef DEBUG write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss", nss write (iout,*) "lenrec_in",lenrec_in @@ -217,19 +217,28 @@ c call flush(iout) if (iret.eq.0) goto 101 call xdrfint_(ixdrf, nss, iret) if (iret.eq.0) goto 101 - do j=1,nss - if (dyn_ss) then - call xdrfint(ixdrf, idssb(j), iret) - call xdrfint(ixdrf, jdssb(j), iret) - idssb(j)=idssb(j)-nres - jdssb(j)=jdssb(j)-nres - else - call xdrfint_(ixdrf, ihpb(j), iret) - if (iret.eq.0) goto 101 - call xdrfint_(ixdrf, jhpb(j), iret) - if (iret.eq.0) goto 101 - endif - enddo + if (dyn_ss) then + do k=1,nss + call xdrfint(ixdrf, idssb(k), iret) + call xdrfint(ixdrf, jdssb(k), iret) + ihpb(k)=iss(idssb(k)-nres)+nres + jhpb(k)=iss(jdssb(k)-nres)+nres +#ifdef DEBUG + write (iout,*) "jj",jj+1," dyn_ss:",idssb(k)-nres, + & jdssb(k)-nres,ihpb(k),jhpb(k) +#endif + enddo + else + do k=1,nss + call xdrfint(ixdrf, ihpb(k), iret) + if (iret.eq.0) goto 101 + call xdrfint(ixdrf, jhpb(k), iret) + if (iret.eq.0) goto 101 +#ifdef DEBUG + write (iout,*) "jj",jj+1," stat_ss:",ihpb(k),jhpb(k) +#endif + enddo + endif call xdrffloat_(ixdrf,reini,iret) if (iret.eq.0) goto 101 call xdrffloat_(ixdrf,refree,iret) @@ -249,17 +258,28 @@ c write (iout,*) "iret",iret c write (iout,*) "nss",nss call flush(iout) if (iret.eq.0) goto 101 - do k=1,nss - if (dyn_ss) then - call xdrfint(ixdrf, idssb(k), iret) - call xdrfint(ixdrf, jdssb(k), iret) + if (dyn_ss) then + do k=1,nss + call xdrfint(ixdrf, idssb(k), iret) + call xdrfint(ixdrf, jdssb(k), iret) + ihpb(k)=iss(idssb(k)-nres)+nres + jhpb(k)=iss(jdssb(k)-nres)+nres +#ifdef DEBUG + write (iout,*) "jj",jj+1," dyn_ss:",idssb(k)-nres, + & jdssb(k)-nres,ihpb(k),jhpb(k) +#endif + enddo else - call xdrfint(ixdrf, ihpb(k), iret) - if (iret.eq.0) goto 101 - call xdrfint(ixdrf, jhpb(k), iret) - if (iret.eq.0) goto 101 + do k=1,nss + call xdrfint(ixdrf, ihpb(k), iret) + if (iret.eq.0) goto 101 + call xdrfint(ixdrf, jhpb(k), iret) + if (iret.eq.0) goto 101 +#ifdef DEBUG + write (iout,*) "jj",jj+1," stat_ss:",ihpb(k),jhpb(k) +#endif + enddo endif - enddo call xdrffloat(ixdrf,reini,iret) if (iret.eq.0) goto 101 call xdrffloat(ixdrf,refree,iret) @@ -391,10 +411,9 @@ c------------------------------------------------------------------------------ include "COMMON.SBRIDGE" include "COMMON.GEO" integer i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib - & nn,nn1,inan,Next,itj,chalen + & nn,nn1,inan,Next,itj double precision etot,energia(0:max_ene) jjj=jjj+1 - chalen=int((nct-nnt+2)/symetr) call int_from_cart1(.false.) do j=nnt+1,nct if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) @@ -402,8 +421,7 @@ c------------------------------------------------------------------------------ if (j.gt.2) then if (itel(j).ne.0 .and. itel(j-1).ne.0) then write (iout,*) "Conformation",jjj,jj+1 - write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j), - & chalen + write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j) write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) @@ -671,17 +689,17 @@ C#define DEBUG write (iout,*) "Reading binary file, record",iii," ii",ii call flush(iout) #endif - if (dyn_ss) then - read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), - & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), +c if (dyn_ss) then +c read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), +c & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), c & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss), - & entfac(ii),rmstb(ii) - else +c & entfac(ii),rmstb(ii) +c else read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss), & entfac(ii),rmstb(ii) - endif +c endif #ifdef DEBUG write (iout,*) ii,iii,ij,entfac(ii) write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres) @@ -738,17 +756,17 @@ c write (iout,*) "Writing binary file, record",iii," ii",ii call flush(iout) #endif - if (dyn_ss) then - write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), - & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), +c if (dyn_ss) then +c write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), +c & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), c & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)) - & entfac(ii),rmstb(ii) - else +c & entfac(ii),rmstb(ii) +c else write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)), & entfac(ii),rmstb(ii) - endif +c endif #ifdef DEBUG write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres) write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres, diff --git a/source/cluster/wham/src-HCD/readrtns.F b/source/cluster/wham/src-HCD/readrtns.F index e9e576f..2cd4ee1 100644 --- a/source/cluster/wham/src-HCD/readrtns.F +++ b/source/cluster/wham/src-HCD/readrtns.F @@ -27,6 +27,7 @@ C read (INP,'(a80)') titel call card_concat(controlcard) + dyn_ss = index(controlcard,"DYN_SS").gt.0 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call readi(controlcard,'TORMODE',tor_mode,0) diff --git a/source/cluster/wham/src-HCD/ssMD.F b/source/cluster/wham/src-HCD/ssMD.F index 9b2908f..a62f7be 100644 --- a/source/cluster/wham/src-HCD/ssMD.F +++ b/source/cluster/wham/src-HCD/ssMD.F @@ -496,8 +496,8 @@ cgrad enddo cgrad enddo do l=1,3 - gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k) - gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k) + gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) enddo return @@ -571,7 +571,7 @@ c Local variables logical found integer i_newnss(1024),displ(0:1024) integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss - + nfgtasks=1 allnss=0 do i=1,ns-1 do j=i+1,ns diff --git a/source/cluster/wham/src-HCD/wrtclust.f b/source/cluster/wham/src-HCD/wrtclust.f index e21494a..0415640 100644 --- a/source/cluster/wham/src-HCD/wrtclust.f +++ b/source/cluster/wham/src-HCD/wrtclust.f @@ -189,7 +189,7 @@ c create InsightII command file for their displaying in different colors c write (iout,*) "cfname",cfname OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED') write (ipdb,'(a,f8.2)') - & "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper + & "REMARK AVERAGE CONFORMATIONS AT TEMPERATURE",temper close (ipdb) I=1 ICON=NCONF(1,1) @@ -227,6 +227,13 @@ c write (iout,*) "ncon_out",ncon_out c(k,ii)=allcart(k,ii,icon) enddo enddo + nss=nss_all(icon) + write (iout,*) "ICON",icon," nss",nss + do k=1,nss + ihpb(k)=ihpb_all(k,icon) + jhpb(k)=jhpb_all(k,icon) + write (iout,*) ihpb(k),jhpb(k) + enddo call center call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel) write (ipdb,'("TER")') @@ -240,6 +247,7 @@ c Average structures and structures closest to average call ave_coord(i) write (ipdb,'(a,i5)') "REMARK CLUSTER",i call center + nss=0 call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel) write (ipdb,'("TER")') if (print_fittest.and.(nsaxs.gt.0 .or. nhpb.gt.0 @@ -530,6 +538,7 @@ c------------------------------------------------------------------------------ include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.INTERACT' + include 'COMMON.SBRIDGE' include 'COMMON.VAR' include 'COMMON.FFIELD' include 'COMMON.TORCNSTR' @@ -570,6 +579,11 @@ c & " Edihcnstr",edihcnstr c(j,i)=allcart(j,i,jconmin) enddo enddo + nss=nss_all(jconmin) + do k=1,nss + ihpb(k)=ihpb_all(k,jconmin) + jhpb(k)=jhpb_all(k,jconmin) + enddo return end c------------------------------------------------------------------------------ @@ -583,6 +597,7 @@ c------------------------------------------------------------------------------ include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.VAR' + include 'COMMON.SBRIDGE' logical non_conv double precision przes(3),obrot(3,3) integer i,ii,j,k,icon,jcon,jconmin,igr,ipermmin @@ -616,6 +631,15 @@ c write (iout,*) "rmsmin",rmsmin," rms",rms c(j,i)=allcart(j,i,jconmin) enddo enddo + nss=nss_all(jconmin) +c write (iout,*) "jconmin",jconmin," nss",nss + call flush(iout) + do k=1,nss + ihpb(k)=ihpb_all(k,jconmin) + jhpb(k)=jhpb_all(k,jconmin) +c write (iout,*) "k",k," ihpb",ihpb(k)," jhpb",jhpb(k) + enddo + call flush(iout) return end c------------------------------------------------------------------------------ -- 1.7.9.5