From 03ac44644b8a2c0bc8f10aeddd4f292c3b2d2e6e Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Tue, 16 Feb 2021 14:54:17 +0100 Subject: [PATCH] Adam's wham update --- source/wham/src-HCD/COMMON.CHAIN | 2 +- source/wham/src-HCD/COMMON.DFA | 2 +- source/wham/src-HCD/COMMON.ENERGIES | 2 +- source/wham/src-HCD/DIMENSIONS | 4 +- source/wham/src-HCD/DIMENSIONS.FREE | 2 +- source/wham/src-HCD/chain_symmetry.F | 3 ++ source/wham/src-HCD/cxread.F | 2 +- source/wham/src-HCD/elecont.f | 63 +++++----------------- source/wham/src-HCD/enecalc1.F | 52 ++++++++++++++++-- source/wham/src-HCD/energy_p_new.F | 5 +- source/wham/src-HCD/geomout.F | 26 ++++----- source/wham/src-HCD/include_unres/COMMON.SBRIDGE | 3 +- source/wham/src-HCD/include_unres/COMMON.SETUP | 4 +- source/wham/src-HCD/initialize_p.F | 2 + source/wham/src-HCD/make_ensemble1.F | 2 +- source/wham/src-HCD/read_constr_homology.F | 4 +- source/wham/src-HCD/secondary.f | 6 +-- source/wham/src-HCD/ssMD.F | 15 ++++-- source/wham/src-HCD/wham_calc1.F | 18 +++---- 19 files changed, 120 insertions(+), 97 deletions(-) diff --git a/source/wham/src-HCD/COMMON.CHAIN b/source/wham/src-HCD/COMMON.CHAIN index 7b79a58..6b453af 100644 --- a/source/wham/src-HCD/COMMON.CHAIN +++ b/source/wham/src-HCD/COMMON.CHAIN @@ -3,7 +3,7 @@ & 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), + 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/wham/src-HCD/COMMON.DFA b/source/wham/src-HCD/COMMON.DFA index 064a7ce..782e8c4 100644 --- a/source/wham/src-HCD/COMMON.DFA +++ b/source/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/wham/src-HCD/COMMON.ENERGIES b/source/wham/src-HCD/COMMON.ENERGIES index 2d40a95..e0ed696 100644 --- a/source/wham/src-HCD/COMMON.ENERGIES +++ b/source/wham/src-HCD/COMMON.ENERGIES @@ -1,4 +1,4 @@ double precision potE(MaxStr_Proc,Max_Parm),entfac(MaxStr_Proc), - & q(MaxQ+2,MaxStr_Proc),enetb(max_ene,MaxStr_Proc,Max_Parm) + & q(MaxQ+6,MaxStr_Proc),enetb(max_ene,MaxStr_Proc,Max_Parm) integer einicheck common /energies/ potE,entfac,q,enetb,einicheck diff --git a/source/wham/src-HCD/DIMENSIONS b/source/wham/src-HCD/DIMENSIONS index fb93081..49e50f5 100644 --- a/source/wham/src-HCD/DIMENSIONS +++ b/source/wham/src-HCD/DIMENSIONS @@ -14,8 +14,8 @@ c parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres c parameter (maxres=250) -c parameter (maxres=1200) - parameter (maxres=10000) + parameter (maxres=1200) +c parameter (maxres=20000) C Max. number of cysteines and other bridging residues integer max_cyst parameter (max_cyst=100) diff --git a/source/wham/src-HCD/DIMENSIONS.FREE b/source/wham/src-HCD/DIMENSIONS.FREE index 7a397d9..c42f2d7 100644 --- a/source/wham/src-HCD/DIMENSIONS.FREE +++ b/source/wham/src-HCD/DIMENSIONS.FREE @@ -3,7 +3,7 @@ integer MaxR,MaxT_h,maxHdim integer MaxSlice parameter (Max_Parm=5) - parameter (MaxQ=4,MaxQ1=MaxQ+2) + parameter (MaxQ=4,MaxQ1=MaxQ+6) parameter(MaxR=8,MaxT_h=36) parameter(MaxSlice=40) parameter(maxHdim=200) diff --git a/source/wham/src-HCD/chain_symmetry.F b/source/wham/src-HCD/chain_symmetry.F index 1406d1d..8c36855 100644 --- a/source/wham/src-HCD/chain_symmetry.F +++ b/source/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/wham/src-HCD/cxread.F b/source/wham/src-HCD/cxread.F index 46867c5..2cfb938 100644 --- a/source/wham/src-HCD/cxread.F +++ b/source/wham/src-HCD/cxread.F @@ -174,7 +174,7 @@ c call flush(iout) c write (iout,*) "Before boxshift" c call flush(iout) c Box shift - call oligomer +c call oligomer c write (iout,*) "After oligomer" c call flush(iout) do i=1,nres diff --git a/source/wham/src-HCD/elecont.f b/source/wham/src-HCD/elecont.f index fb105a4..86db2df 100644 --- a/source/wham/src-HCD/elecont.f +++ b/source/wham/src-HCD/elecont.f @@ -18,10 +18,11 @@ & eesij,ees,evdw,ene, rij,zj_temp,xj_temp,yj_temp, & sscale,sscagrad,dist_temp,xj_safe,yj_safe,zj_safe,dist_init double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2), - & appc(2,2),bppc(2,2) + & appc(2,2),bppc(2,2),epp_(2,2),rpp_(2,2) double precision elcutoff,elecutoff_14 integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap double precision econt(maxcont) + double precision boxshift * * Load the constants of peptide bond - peptide bond interactions. * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g. @@ -29,8 +30,8 @@ * * as of 7/06/91. * -c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ -c data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ +c data epp_ / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ + data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ data elpp6c /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/ data elpp3c / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/ data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/ @@ -40,7 +41,7 @@ c data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ & "Constants of electrostatic interaction energy expression." do i=1,2 do j=1,2 - rri=rpp(i,j)**6 + rri=rpp_(i,j)**6 appc(i,j)=epp(i,j)*rri*rri bppc(i,j)=-2.0*epp(i,j)*rri ael6c(i,j)=elpp6c(i,j)*4.2**6 @@ -62,12 +63,8 @@ c data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ xmedi=xi+0.5*dxi ymedi=yi+0.5*dyi zmedi=zi+0.5*dzi - xmedi=mod(xmedi,boxxsize) - if (xmedi.lt.0) xmedi=xmedi+boxxsize - ymedi=mod(ymedi,boxysize) - if (ymedi.lt.0) ymedi=ymedi+boxysize - zmedi=mod(zmedi,boxzsize) - if (zmedi.lt.0) zmedi=zmedi+boxzsize + call to_box(xmedi,ymedi,zmedi) +c write (iout,*) "i",xmedi,ymedi,zmedi do 4 j=i+2,ien-1 jj=iperm(j,ipermmin) ind=ind+1 @@ -86,46 +83,13 @@ c data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ xj=c(1,jj)+0.5*dxj yj=c(2,jj)+0.5*dyj zj=c(3,jj)+0.5*dzj - xj=mod(xj,boxxsize) - if (xj.lt.0) xj=xj+boxxsize - yj=mod(yj,boxysize) - if (yj.lt.0) yj=yj+boxysize - zj=mod(zj,boxzsize) - if (zj.lt.0) zj=zj+boxzsize - dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - isubchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 - xj=xj_safe+xshift*boxxsize - yj=yj_safe+yshift*boxysize - zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - isubchap=1 - endif - enddo - enddo - enddo - if (isubchap.eq.1) then - xj=xj_temp-xmedi - yj=yj_temp-ymedi - zj=zj_temp-zmedi - else - xj=xj_safe-xmedi - yj=yj_safe-ymedi - zj=zj_safe-zmedi - endif +c write (iout,*) "j",xj,yj,zj + call to_box(xj,yj,zj) + xj=boxshift(xj-xmedi,boxxsize) + yj=boxshift(yj-ymedi,boxysize) + zj=boxshift(zj-zmedi,boxzsize) +c write (iout,*) "j",xj,yj,zj rij=xj*xj+yj*yj+zj*zj - sss=sscale(sqrt(rij)) - sssgrad=sscagrad(sqrt(rij)) rrmij=1.0/(xj*xj+yj*yj+zj*zj) rmij=sqrt(rrmij) r3ij=rrmij*rmij @@ -152,6 +116,7 @@ c data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ endif ees=ees+eesij evdw=evdw+evdwij*sss +c write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij 4 continue 1 continue if (lprint) then diff --git a/source/wham/src-HCD/enecalc1.F b/source/wham/src-HCD/enecalc1.F index b64de48..2edf349 100644 --- a/source/wham/src-HCD/enecalc1.F +++ b/source/wham/src-HCD/enecalc1.F @@ -46,13 +46,31 @@ c double precision tole /1.0d-1/ double precision tt integer snk_p(MaxR,MaxT_h,Max_parm) logical lerr + integer ncont,icont(2,maxcont),isecstr(maxres) character*256 bprotfile_temp + double precision totlength call opentmp(islice,ientout,bprotfile_temp) iii=0 ii=0 errmsg_count=0 c write (iout,*) "enecalc: nparmset ",nparmset c write (iout,*) "enecalc: tormode ",tor_mode + write (iout,*) "ns",ns," dyn_ss",dyn_ss,(iss(i),i=1,ns) + if (ns.gt.0.and.dyn_ss) then + do i=nss+1,nhpb + ihpb(i-nss)=ihpb(i) + jhpb(i-nss)=jhpb(i) + forcon(i-nss)=forcon(i) + dhpb(i-nss)=dhpb(i) + enddo + nhpb=nhpb-nss + nss=0 + call hpb_partition + do i=1,ns + dyn_ss_mask(iss(i))=.true. + enddo + endif + write (iout,*) "dyn_ss_mask",(dyn_ss_mask(i),i=1,nres) #ifdef MPI do iparm=1,nParmSet do ib=1,nT_h(iparm) @@ -94,7 +112,7 @@ c write (iout,*) "enecalc: tormode ",tor_mode anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3) q(nQ+1,iii+1)=rmsnat(iii+1,ipermin) endif - q(nQ+2,iii+1)=gyrate(iii+1) +c write (iout,*) iii+1,q(nQ+3,iii+1),q(nQ+4,iii+1),q(nQ+5,iii+1) c fT=T0*beta_h(ib,ipar)*1.987D-3 c ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)) if (rescale_mode.eq.1) then @@ -158,9 +176,9 @@ c & " kfac",kfac,"quot",quot," fT",fT & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, & wtor_d,wsccor,wbond #endif -C write (iout,*) "tuz przed energia" +c write (iout,*) "tuz przed energia" call etotal(energia(0),fT) -C write (iout,*) "tuz za energia" +c write (iout,*) "tuz za energia" #ifdef DEBUG write (iout,*) "Conformation",i write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres), @@ -215,8 +233,8 @@ c & eini-energia(0) write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres), & ((c(l,k+nres),l=1,3),k=nnt,nct) c call intout - call pdbout(indstart(me1)+iii, - & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0) +c call pdbout(indstart(me1)+iii, +c & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0) call enerprint(energia(0),fT) errmsg_count=errmsg_count+1 if (errmsg_count.gt.maxerrmsg_count) @@ -260,6 +278,30 @@ c call enerprint(energia(0),fT) endif enddo ! iparm + q(nQ+2,iii+1)=gyrate(iii+1) +c 8/28/2020 Adam - determine the fraction of secondary structures. + call elecont(.false.,ncont,icont,nnt,nct-1,1) + call secondary2(.false.,.false.,ncont,icont,isecstr) +#ifdef DEBUG + write (iout,*) "secondary structure" + write (iout,'(80i1)') (isecstr(k),k=1,nres) +#endif + q(nQ+3,iii+1)=0.0d0 + q(nQ+4,iii+1)=0.0d0 + q(nQ+5,iii+1)=0.0d0 + totlength=0.0d0 + do k=nnt,nct + if (itype(k).eq.ntyp1) cycle + totlength=totlength+1.0d0 + l=isecstr(k) + q(nQ+3+l,iii+1)=q(nQ+3+l,iii+1)+1.0d0 + enddo + q(nQ+3,iii+1)=q(nQ+3,iii+1)/totlength + q(nQ+4,iii+1)=q(nQ+4,iii+1)/totlength + q(nQ+5,iii+1)=q(nQ+5,iii+1)/totlength +c write (iout,*) "iii",iii," nssbond",nssbond,nss +c q(nQ+6,iii+1)=nssbond + q(nQ+6,iii+1)=nss iii=iii+1 if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) diff --git a/source/wham/src-HCD/energy_p_new.F b/source/wham/src-HCD/energy_p_new.F index e72d558..3a83918 100644 --- a/source/wham/src-HCD/energy_p_new.F +++ b/source/wham/src-HCD/energy_p_new.F @@ -56,7 +56,7 @@ C call set_shield_fac2 endif call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) -C write(iout,*) 'po eelec' +c write(iout,*) 'po eelec eello_turn4',eello_turn4 C Calculate excluded-volume interaction energy between peptide groups C and side chains. @@ -2273,6 +2273,7 @@ c write(iout,*) "JESTEM W PETLI" call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) +c write (iout,*) "i",i," eello_turn4",eello_turn4 #ifdef FOURBODY num_cont_hb(i)=num_conti #endif @@ -3581,6 +3582,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/wham/src-HCD/geomout.F b/source/wham/src-HCD/geomout.F index 097040f..ed33cc7 100644 --- a/source/wham/src-HCD/geomout.F +++ b/source/wham/src-HCD/geomout.F @@ -9,7 +9,9 @@ include 'COMMON.HEADER' include 'COMMON.SBRIDGE' character*50 tytul - character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/ + 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'/ dimension ica(maxres) write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)') & ii,temp,rmsdev @@ -20,24 +22,28 @@ iatom=0 ichain=1 ires=0 + iti_prev=0 do i=nnt,nct iti=itype(i) if (iti.eq.ntyp1) then - ichain=ichain+1 ires=0 - write (ipdb,'(a)') 'TER' + if (iti_prev.ne.ntyp1) then + write (ipdb,'(a)') 'TER' + ichain=ichain+1 + endif else ires=ires+1 iatom=iatom+1 ica(i)=iatom write (ipdb,10) iatom,restyp(iti),chainid(ichain), - & ires,(c(j,i),j=1,3) + & ires,(c(j,i),j=1,3),1.0d0 if (iti.ne.10) then iatom=iatom+1 write (ipdb,20) iatom,restyp(iti),chainid(ichain), - & ires,(c(j,nres+i),j=1,3) + & ires,(c(j,nres+i),j=1,3),1.0d0 endif endif + iti_prev=iti enddo write (ipdb,'(a)') 'TER' do i=nnt,nct-1 @@ -54,15 +60,11 @@ write (ipdb,30) ica(nct),ica(nct)+1 endif do i=1,nss - if (dyn_ss) then - write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1 - else - write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1 - endif + write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1 enddo write (ipdb,'(a6)') 'ENDMDL' - 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3) - 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3) + 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2) + 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2) 30 FORMAT ('CONECT',8I5) return end diff --git a/source/wham/src-HCD/include_unres/COMMON.SBRIDGE b/source/wham/src-HCD/include_unres/COMMON.SBRIDGE index a313d8f..d2a41e1 100644 --- a/source/wham/src-HCD/include_unres/COMMON.SBRIDGE +++ b/source/wham/src-HCD/include_unres/COMMON.SBRIDGE @@ -24,8 +24,9 @@ & link_end_peak double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss logical dyn_ss,dyn_ss_mask + integer nssbond common /dyn_ssbond/ dtriss,atriss,btriss,ctriss,Ht, & dyn_ssbond_ij(max_cyst,max_cyst), - & idssb(maxss),jdssb(maxss) + & idssb(maxss),jdssb(maxss),nssbond common /dyn_ss_logic/ & dyn_ss,dyn_ss_mask(maxres) diff --git a/source/wham/src-HCD/include_unres/COMMON.SETUP b/source/wham/src-HCD/include_unres/COMMON.SETUP index 5039116..1edc7c3 100644 --- a/source/wham/src-HCD/include_unres/COMMON.SETUP +++ b/source/wham/src-HCD/include_unres/COMMON.SETUP @@ -1,14 +1,14 @@ integer king,idint,idreal,idchar,is_done parameter (king=0,idint=1105,idreal=1729,idchar=1597,is_done=1) integer me,cg_rank,fg_rank,fg_rank1,nodes,Nprocs,nfgtasks,kolor, - & koniec(0:maxprocs-1),WhatsUp,ifinish(maxprocs-1),CG_COMM,FG_COMM, + & koniec(0:maxprocs-1),ifinish(maxprocs-1),CG_COMM,FG_COMM, & FG_COMM1,CONT_FROM_COMM,CONT_TO_COMM,lentyp(0:maxprocs-1), & kolor1,key1,nfgtasks1,MyRank, & max_gs_size logical yourjob, finished, cgdone common/setup/me,MyRank,cg_rank,fg_rank,fg_rank1,nodes,Nprocs, & nfgtasks,nfgtasks1, - & max_gs_size,kolor,koniec,WhatsUp,ifinish,CG_COMM,FG_COMM, + & max_gs_size,kolor,koniec,ifinish,CG_COMM,FG_COMM, & FG_COMM1,CONT_FROM_COMM,CONT_TO_COMM,lentyp integer MPI_UYZ,MPI_UYZGRAD,MPI_MU,MPI_MAT1,MPI_MAT2, & MPI_THET,MPI_GAM, diff --git a/source/wham/src-HCD/initialize_p.F b/source/wham/src-HCD/initialize_p.F index 8217663..0cf093a 100644 --- a/source/wham/src-HCD/initialize_p.F +++ b/source/wham/src-HCD/initialize_p.F @@ -22,6 +22,7 @@ C include "COMMON.NAMES" include "COMMON.TIME1" include "COMMON.TORCNSTR" + include "COMMON.SETUP" C C The following is just to define auxiliary variables used in angle conversion C @@ -237,6 +238,7 @@ C Set timers and counters for the respective routines n_viol = 0 n_gviol = 0 n_map = 0 + nfgtasks = 1 #ifndef SPLITELE nprint_ene=nprint_ene-1 #endif diff --git a/source/wham/src-HCD/make_ensemble1.F b/source/wham/src-HCD/make_ensemble1.F index a07dbeb..f18dd2b 100644 --- a/source/wham/src-HCD/make_ensemble1.F +++ b/source/wham/src-HCD/make_ensemble1.F @@ -371,7 +371,7 @@ c write (iout,*) "qfree",qfree & ctemper(:ilen(ctemper))//"pdb" endif open(ipdb,file=pdbname) - write (iout,*) "Before reading nlist",nlist +c write (iout,*) "Before reading nlist",nlist do i=1,nlist read (ientout,rec=iperm(i)) & ((csingle(l,k),l=1,3),k=1,nres), diff --git a/source/wham/src-HCD/read_constr_homology.F b/source/wham/src-HCD/read_constr_homology.F index fa81b80..03d7968 100644 --- a/source/wham/src-HCD/read_constr_homology.F +++ b/source/wham/src-HCD/read_constr_homology.F @@ -259,8 +259,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/wham/src-HCD/secondary.f b/source/wham/src-HCD/secondary.f index 4088831..d2e9271 100644 --- a/source/wham/src-HCD/secondary.f +++ b/source/wham/src-HCD/secondary.f @@ -401,8 +401,8 @@ cd write (iout,*) '------- looking for parallel beta -----------' do i=1,ncont i1=icont(1,i) j1=icont(2,i) - if (i1.ge.nstart_sup .and. i1.le.nend_sup - & .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then +c if (i1.ge.nstart_sup .and. i1.le.nend_sup +c & .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then cd write (iout,*) "parallel",i1,j1 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then ii1=i1 @@ -469,7 +469,7 @@ cd write (iout,*) i1,j1,not_done & "SetNeigh",ii1-1,i1-1,jj1-1,j1-1 endif endif - endif +c endif endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup enddo diff --git a/source/wham/src-HCD/ssMD.F b/source/wham/src-HCD/ssMD.F index 4ce1b3d..d9b9df7 100644 --- a/source/wham/src-HCD/ssMD.F +++ b/source/wham/src-HCD/ssMD.F @@ -144,13 +144,13 @@ c-------TESTING CODE double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi c-------END TESTING CODE - + nssbond=0 i=resi j=resj ici=icys(i) icj=icys(j) if (ici.eq.0 .or. icj.eq.0) then - write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)') + write (iout,'(a,i5,2a,a3,i5,5h and ,a3,i5)') & "Attempt to create", & " a disulfide link between non-cysteine residues ",restyp(i),i, & restyp(j),j @@ -276,6 +276,8 @@ c Stop and plot energy and derivative as a function of distance & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps c-------END TESTING CODE +c write (iout,'(2(a,i5),4(a,f7.2))') "resi",resi," resj",resj, +c & " ljxm",ljxm," ljxs",ljxs," ssxm",ssxm," rij",rij if (rij.gt.ljxm) then havebond=.false. ljd=rij-ljXs @@ -298,6 +300,8 @@ c-------END TESTING CODE & -2.0D0*alf12*eps3der+sigder*sigsq_om12 else if (rij.lt.ssxm) then havebond=.true. + nssbond=nssbond+1 +c write (iout,*) "ssMD: nssbond",nssbond ssd=rij-ssXs eij=ssA*ssd*ssd+ssB*ssd+ssC eij=eij*sss @@ -309,6 +313,7 @@ c-------END TESTING CODE eom2= 2*akth*deltat2+pom1-om1*pom2 eom12=pom2 else +c nssbond=nssbond+1 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi d_ssxm(1)=0.5D0*akct/ssA @@ -497,8 +502,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 @@ -556,7 +561,6 @@ c Includes include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' -C include 'COMMON.SETUP' #ifndef CLUST #ifndef WHAM C include 'COMMON.MD' @@ -572,6 +576,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 diff --git a/source/wham/src-HCD/wham_calc1.F b/source/wham/src-HCD/wham_calc1.F index 7e4512d..2d1d661 100644 --- a/source/wham/src-HCD/wham_calc1.F +++ b/source/wham/src-HCD/wham_calc1.F @@ -639,7 +639,7 @@ c & " WHAM_COMM",WHAM_COMM sumE_p(i,iparm)=0.0d0 sumEbis_p(i,iparm)=0.0d0 sumEsq_p(i,iparm)=0.0d0 - do j=1,nQ+2 + do j=1,nQ+6 sumQ_p(j,i,iparm)=0.0d0 sumQsq_p(j,i,iparm)=0.0d0 sumEQ_p(j,i,iparm)=0.0d0 @@ -654,7 +654,7 @@ c & " WHAM_COMM",WHAM_COMM sumE(i,iparm)=0.0d0 sumEbis(i,iparm)=0.0d0 sumEsq(i,iparm)=0.0d0 - do j=1,nQ+2 + do j=1,nQ+6 sumQ(j,i,iparm)=0.0d0 sumQsq(j,i,iparm)=0.0d0 sumEQ(j,i,iparm)=0.0d0 @@ -854,7 +854,7 @@ c call restore_parm(iparm) sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight - do j=1,nQ+2 + do j=1,nQ+6 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) @@ -865,7 +865,7 @@ c call restore_parm(iparm) sumE(k,iparm)=sumE(k,iparm)+etot*weight sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight - do j=1,nQ+2 + do j=1,nQ+6 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight sumEQ(j,k,iparm)=sumEQ(j,k,iparm) @@ -1085,7 +1085,7 @@ c call restore_parm(iparm) & sumW(i,iparm) sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2) - do j=1,nQ+2 + do j=1,nQ+6 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm) sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) & -sumQ(j,i,iparm)**2 @@ -1096,15 +1096,15 @@ c call restore_parm(iparm) & (startGridT+i*delta_T))+potEmin write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T, & sumW(i,iparm),sumE(i,iparm) - write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2) + write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+6) write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm), - & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2) + & (sumQsq(j,i,iparm),j=1,nQ+6),(sumEQ(j,i,iparm),j=1,nQ+6) write (iout,*) write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T, & sumW(i,iparm),sumE(i,iparm) - write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2) + write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+6) write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm), - & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2) + & (sumQsq(j,i,iparm),j=1,nQ+6),(sumEQ(j,i,iparm),j=1,nQ+6) write (34,*) call flush(34) enddo -- 1.7.9.5