From ae315105cef83bbcab70e2778ef92459690ee784 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Sun, 19 Jul 2020 19:32:07 +0200 Subject: [PATCH] changes --- source/unres/src-HCD-5D/COMMON.DFA | 2 +- source/unres/src-HCD-5D/DIMENSIONS | 7 ++-- source/unres/src-HCD-5D/boxshift.f | 12 ++++++- source/unres/src-HCD-5D/chain_symmetry.F | 3 ++ source/unres/src-HCD-5D/dfa.F | 1 + source/unres/src-HCD-5D/energy_p_new_barrier.F | 3 +- source/unres/src-HCD-5D/readpdb-mult.F | 43 +++++++----------------- source/unres/src-HCD-5D/readrtns_CSA.F | 18 +++++++--- source/unres/src-HCD-5D/ssMD.F | 29 +++++++++++++--- source/wham/src-HCD/bxread.F | 2 +- source/wham/src-HCD/cxread.F | 2 +- source/wham/src-HCD/enecalc1.F | 4 +-- source/wham/src-HCD/readrtns.F | 4 +-- source/wham/src-HCD/xread.F | 2 +- 14 files changed, 80 insertions(+), 52 deletions(-) diff --git a/source/unres/src-HCD-5D/COMMON.DFA b/source/unres/src-HCD-5D/COMMON.DFA index 818c024..51a7af7 100644 --- a/source/unres/src-HCD-5D/COMMON.DFA +++ b/source/unres/src-HCD-5D/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/unres/src-HCD-5D/DIMENSIONS b/source/unres/src-HCD-5D/DIMENSIONS index ed21dfe..4236776 100644 --- a/source/unres/src-HCD-5D/DIMENSIONS +++ b/source/unres/src-HCD-5D/DIMENSIONS @@ -16,7 +16,7 @@ C Max. number of coarse-grain processors parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres - parameter (maxres=10000) + parameter (maxres=20000) C Max. number of AA residues per chain integer maxres_chain parameter (maxres_chain=1200) @@ -53,14 +53,15 @@ C Max. number of contacts per residue c parameter (maxconts=50) C Max. number of interactions within cutoff per residue integer maxint_res - parameter (maxint_res=200) + parameter (maxint_res=250) C Max. number od residues within distance cufoff from a given residue to C include in template-based/contact distance restraints. integer maxcont_res parameter (maxcont_res=200) C Max. number of distance/contact-distance restraints integer maxdim_cont - parameter (maxdim_cont=maxres*maxcont_res) +c parameter (maxdim_cont=maxres*maxcont_res) + parameter (maxdim_cont=maxres*1000) C Number of AA types (at present only natural AA's will be handled integer ntyp,ntyp1 parameter (ntyp=24,ntyp1=ntyp+1) diff --git a/source/unres/src-HCD-5D/boxshift.f b/source/unres/src-HCD-5D/boxshift.f index 29d3406..070d9c9 100644 --- a/source/unres/src-HCD-5D/boxshift.f +++ b/source/unres/src-HCD-5D/boxshift.f @@ -72,11 +72,17 @@ c-------------------------------------------------------------------------- subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi) implicit none include 'DIMENSIONS' + include 'COMMON.IOUNITS' include 'COMMON.CHAIN' double precision xi,yi,zi,sslipi,ssgradlipi double precision fracinbuf double precision sscalelip,sscagradlip - +#define DEBUG +#ifdef DEBUG + write (iout,*) "bordlipbot",bordlipbot," bordliptop",bordliptop + write (iout,*) "buflipbot",buflipbot," lipbufthick",lipbufthick + write (iout,*) "xi yi zi",xi,yi,zi +#endif if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then C the energy transfer exist if (zi.lt.buflipbot) then @@ -97,5 +103,9 @@ C lipbufthick is thickenes of lipid buffore sslipi=0.0d0 ssgradlipi=0.0 endif +#ifdef DEBUG + write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi +#endif +#undef DEBUG return end diff --git a/source/unres/src-HCD-5D/chain_symmetry.F b/source/unres/src-HCD-5D/chain_symmetry.F index 1406d1d..8c36855 100644 --- a/source/unres/src-HCD-5D/chain_symmetry.F +++ b/source/unres/src-HCD-5D/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/unres/src-HCD-5D/dfa.F b/source/unres/src-HCD-5D/dfa.F index 3982b6d..1af0b44 100644 --- a/source/unres/src-HCD-5D/dfa.F +++ b/source/unres/src-HCD-5D/dfa.F @@ -368,6 +368,7 @@ C END OF BETA RESTRAINT edfadis=0 gdfad=0.0d0 +c write (2,*) "edfad",idfadis_start,idfadis_end do i=idfadis_start,idfadis_end iatm1=idislis(1,i)+ishiftca diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F index 6d5b25f..6d6a817 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -2044,7 +2044,8 @@ c write(iout,*) "PO ZWYKLE", evdwij if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & 'evdw',i,j,evdwij,' ss' C triple bond artifac removal - do k=j+1,iend(i,iint) +c do k=j+1,iend(i,iint) + do k=j+1,nct C search over all next residues if (dyn_ss_mask(k)) then C check if they are cysteins diff --git a/source/unres/src-HCD-5D/readpdb-mult.F b/source/unres/src-HCD-5D/readpdb-mult.F index 76cb6b6..41fe7f6 100644 --- a/source/unres/src-HCD-5D/readpdb-mult.F +++ b/source/unres/src-HCD-5D/readpdb-mult.F @@ -64,9 +64,9 @@ c call flush(iout) iterter(ires_old-1)=1 itype(ires_old)=ntyp1 iterter(ires_old)=1 -c ishift1=ishift1+1 + ishift1=ishift1+1 ibeg=2 - write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg + write (iout,*) "Chain ended",ires,ishift,ires_old if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -95,8 +95,8 @@ c write (iout,*) "! ",atom," !",ires read (card(18:20),'(a3)') res c write (iout,*) "ires",ires,ires-ishift+ishift1, c & " ires_old",ires_old -c write (iout,*) "ishift",ishift," ishift1",ishift1 -c write (iout,*) "IRES",ires-ishift+ishift1,ires_old +c write (iout,*) "ishift",ishift," ishift1",ishift1 +c write (iout,*) "IRES",ires-ishift+ishift1,ires_old if (ires-ishift+ishift1.ne.ires_old) then ! Calculate the CM of the preceding residue. ! if (ibeg.eq.0) call sccenter(ires,iii,sccor) @@ -115,7 +115,6 @@ c write (iout,'(i5,3f10.5)') ires,(sccor(j,iii),j=1,3) sccalc=.true. endif ! Start new residue. -c write (iout,*) "ibeg",ibeg if (res.eq.'Cl-' .or. res.eq.'Na+') then ires=ires_old cycle @@ -134,13 +133,10 @@ c write (iout,*) "BEG ires",ires else if (ibeg.eq.2) then ! Start a new chain ishift=-ires_old+ires-1 !!!!! -c ishift1=ishift1-1 !!!!! -c write (iout,*) "New chain started",ires,ires_old,ishift, -c & ishift1 + ishift1=ishift1-1 !!!!! +c write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 - write (iout,*) "New chain started ires",ires ires_old=ires -c ires=ires_old+1 ibeg=0 else ishift=ishift-(ires-ishift+ishift1-ires_old-1) @@ -163,8 +159,7 @@ c write (2,*) "ires",ires," res ",res!," ity"!,ity if (atom.eq.'CA' .or. atom.eq.'CH3' .or. & res.eq.'NHE'.and.atom(:2).eq.'HN') then read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) -c write (iout,*) "backbone ",atom -c write (iout,*) ires,res,(c(j,ires),j=1,3) +! write (iout,*) "backbone ",atom #ifdef DEBUG write (iout,'(i6,i3,2x,a,3f8.3)') & ires,itype(ires),res,(c(j,ires),j=1,3) @@ -234,7 +229,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + c(j,i)=c(j,i+1)-1.9d0*e2(j) enddo else !unres_pdb do j=1,3 @@ -417,9 +412,8 @@ C and convert the peptide geometry into virtual-chain geometry. character*3 seq,res character*5 atom character*80 card - double precision sccor(3,50) + double precision sccor(3,20) integer rescode,iterter(maxres) - logical zero do i=1,maxres iterter(i)=0 enddo @@ -587,7 +581,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + c(j,i)=c(j,i+1)-1.9d0*e2(j) enddo else !unres_pdb do j=1,3 @@ -622,7 +616,7 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + c(j,nres)=c(j,nres-1)-1.9d0*e2(j) enddo else do j=1,3 @@ -654,7 +648,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0) + c(j,1)=c(j,2)-1.9d0*e2(j) enddo else do j=1,3 @@ -682,18 +676,6 @@ C Calculate internal coordinates. & (c(j,ires+nres),j=1,3) enddo endif - zero=.false. - do ires=1,nres - zero=zero.or.itype(ires).eq.0 - enddo - if (zero) then - write (iout,'(2a)') "Gaps in PDB coordinates detected;", - & " look for ZERO in the control output above." - write (iout,'(2a)') "Repair the PDB file using MODELLER", - & " or other softwared and resubmit." - call flush(iout) - stop - endif C Calculate internal coordinates. call int_from_cart(.true.,out_template_coord) call sc_loc_geom(.false.) @@ -701,7 +683,6 @@ C Calculate internal coordinates. thetaref(i)=theta(i) phiref(i)=phi(i) enddo - dc(:,0)=c(:,1) do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index eeaf74c..4fbc0f1 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -741,7 +741,7 @@ C integer ilen external ilen integer iperm,tperm - integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp + integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp,itemp double precision sumv C C Read PDB structure if applicable @@ -1184,8 +1184,7 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb 335 continue unres_pdb=.false. nres_temp=nres -c call readpdb - call readpdb_template(nmodel_start+1) + call readpdb close(ipdbin) if (nres.ge.nres_temp) then nmodel_start=nmodel_start+1 @@ -1196,11 +1195,22 @@ c call readpdb enddo enddo else - if (me.eq.king .or. .not. out1file) +c itemp=nres +c nres=nres_temp +c call gen_rand_conf(itemp,*115) +c nmodel_start=nmodel_start+1 +c do i=1,2*nres +c do j=1,3 +c chomo(j,i,nmodel_start)=c(j,i) +c enddo +c enddo +c goto 116 + 115 if (me.eq.king .or. .not. out1file) & write (iout,'(a,2i5,1x,a)') & "Different number of residues",nres_temp,nres, & " model skipped." endif + 116 continue nres=nres_temp enddo 332 continue diff --git a/source/unres/src-HCD-5D/ssMD.F b/source/unres/src-HCD-5D/ssMD.F index 26807a0..67e9f10 100644 --- a/source/unres/src-HCD-5D/ssMD.F +++ b/source/unres/src-HCD-5D/ssMD.F @@ -100,6 +100,7 @@ C----------------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.NAMES' + include 'COMMON.SPLITELE' #ifndef CLUST #ifndef WHAM include 'COMMON.MD' @@ -153,6 +154,8 @@ c-------END TESTING CODE j=resj ici=icys(i) icj=icys(j) +c write (iout,*) "dyn_ssbond",resi,resj,ici,icj +c call flush(iout) if (ici.eq.0 .or. icj.eq.0) then #ifdef MPI write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)') @@ -177,6 +180,8 @@ c-------END TESTING CODE yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) +c write (iout,*) "After to_box i",xi,yi,zi +c call flush(iout) C define scaling factor for lipids C if (positi.le.0) positi=positi+boxzsize @@ -184,12 +189,18 @@ C print *,i C first for peptide groups c for each residue check if it is in lipid or lipid water border area call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) +c write (iout,*) "After lipid_layer" +c call flush(iout) itypj=itype(j) xj=c(1,nres+j) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) +c write (iout,*) "After to_box j",xj,yj,zj +c call flush(iout) call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) +c write (iout,*) "After lipid_layer" +c call flush(iout) aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 @@ -197,6 +208,8 @@ c for each residue check if it is in lipid or lipid water border area xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) +c write (iout,*) "After boxshift" +c call flush(iout) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -214,8 +227,8 @@ c for each residue check if it is in lipid or lipid water border area rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse - sss=sscale((1.0d0/rij)/sigma(itypi,itypj)) - sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) + sss=sscale((1.0d0/rij)/sigma(itypi,itypj),r_cut_int) + sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj),r_cut_int) c The following are set in sc_angular c erij(1)=xj*rij c erij(2)=yj*rij @@ -223,7 +236,11 @@ c erij(3)=zj*rij c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3) c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3) c om12=dxi*dxj+dyi*dyj+dzi*dzj +c write (iout,*) "Calling sc_angular" +c call flush(iout) call sc_angular +c write (iout,*) "After sc_angular" +c call flush(iout) rij=1.0D0/rij ! Reset this so it makes sense sig0ij=sigma(itypi,itypj) @@ -265,6 +282,8 @@ c-------END TESTING CODE c-------TESTING CODE c Stop and plot energy and derivative as a function of distance +c write (iout,*) "checkstop",checkstop +c call flush(iout) if (checkstop) then ssm=ssC-0.25D0*ssB*ssB/ssA ljm=-0.25D0*ljB*bb/aa @@ -447,6 +466,8 @@ c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE endif +c write (iout,*) "havebond",havebond +c call flush(iout) if (havebond) then #ifndef CLUST #ifndef WHAM @@ -509,8 +530,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 diff --git a/source/wham/src-HCD/bxread.F b/source/wham/src-HCD/bxread.F index c459499..a85bff4 100644 --- a/source/wham/src-HCD/bxread.F +++ b/source/wham/src-HCD/bxread.F @@ -20,7 +20,7 @@ include "COMMON.FREE" include "COMMON.SBRIDGE" real*4 csingle(3,maxres2) - character*64 nazwa,bprotfile_temp + character*256 nazwa,bprotfile_temp character*3 liczba integer i,is,ie,j,ii,jj,k,kk,l,ll,mm,if integer nrec,nlines,iscor,islice diff --git a/source/wham/src-HCD/cxread.F b/source/wham/src-HCD/cxread.F index 36ef6e6..46867c5 100644 --- a/source/wham/src-HCD/cxread.F +++ b/source/wham/src-HCD/cxread.F @@ -17,7 +17,7 @@ include 'COMMON.VAR' include 'COMMON.GEO' include 'COMMON.PROT' - character*64 nazwa,bprotfile_temp + character*256 nazwa,bprotfile_temp real*4 rtime,rpotE,ruconst,rt_bath,rprop(maxQ) double precision time integer iret,itmp,itraj,ntraj diff --git a/source/wham/src-HCD/enecalc1.F b/source/wham/src-HCD/enecalc1.F index 60addc7..b64de48 100644 --- a/source/wham/src-HCD/enecalc1.F +++ b/source/wham/src-HCD/enecalc1.F @@ -46,7 +46,7 @@ c double precision tole /1.0d-1/ double precision tt integer snk_p(MaxR,MaxT_h,Max_parm) logical lerr - character*64 bprotfile_temp + character*256 bprotfile_temp call opentmp(islice,ientout,bprotfile_temp) iii=0 ii=0 @@ -368,7 +368,7 @@ c------------------------------------------------------------------------------ include "COMMON.CONTACTS1" character*64 nazwa character*80 bxname,cxname - character*64 bprotfile_temp + character*256 bprotfile_temp character*3 liczba,licz character*2 licz2 integer i,itj,ii,iii,j,k,l diff --git a/source/wham/src-HCD/readrtns.F b/source/wham/src-HCD/readrtns.F index bca3771..c733c37 100644 --- a/source/wham/src-HCD/readrtns.F +++ b/source/wham/src-HCD/readrtns.F @@ -412,7 +412,7 @@ c------------------------------------------------------------------------------- include "COMMON.PROTFILES" include "COMMON.PROT" include "COMMON.FREE" - character*64 bprotfile_temp + character*256 bprotfile_temp character*3 liczba,liczba2 character*2 liczba1 integer iunit,islice @@ -472,7 +472,7 @@ c------------------------------------------------------------------------------- include "COMMON.SBRIDGE" include "COMMON.OBCINKA" real*4 csingle(3,maxres2) - character*64 nazwa,bprotfile_temp + character*256 nazwa,bprotfile_temp character*3 liczba character*2 liczba1 integer i,j,ii,jj(maxslice),k,kk(maxslice),l, diff --git a/source/wham/src-HCD/xread.F b/source/wham/src-HCD/xread.F index ac35de1..b397a6f 100644 --- a/source/wham/src-HCD/xread.F +++ b/source/wham/src-HCD/xread.F @@ -23,7 +23,7 @@ include "COMMON.SBRIDGE" include "COMMON.OBCINKA" real*4 csingle(3,maxres2) - character*64 nazwa,bprotfile_temp + character*256 nazwa,bprotfile_temp integer i,j,k,l,ii,jj(maxslice),kk(maxslice),ll(maxslice), & mm(maxslice) integer iscor,islice,islice1,slice -- 1.7.9.5