From a1081b395546875bbbdb59cdd25f52d974adf2a6 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Sun, 24 May 2015 14:46:48 +0200 Subject: [PATCH] Corrected WHAM and UNRES reading in SS dyn... --- source/unres/src_MD-M/parmread.F | 2 +- source/wham/src-M/geomout.F | 10 +-- source/wham/src-M/include_unres/COMMON.SBRIDGE | 4 +- source/wham/src-M/parmread.F | 91 ++++++++++++++++++------ source/wham/src-M/ssMD.F | 7 +- 5 files changed, 81 insertions(+), 33 deletions(-) diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 3296e45..c841df2 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -1256,7 +1256,7 @@ c lprint=.false. C C Define the constants of the disulfide bridge C - ebr=-12.00D0 +C ebr=-12.00D0 c c Old arbitrary potential - commented out. c diff --git a/source/wham/src-M/geomout.F b/source/wham/src-M/geomout.F index c9110ae..6747e76 100644 --- a/source/wham/src-M/geomout.F +++ b/source/wham/src-M/geomout.F @@ -54,11 +54,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 +C if (dyn_ss) then +C write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1 +C else +C write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1 +C endif enddo write (ipdb,'(a6)') 'ENDMDL' 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3) diff --git a/source/wham/src-M/include_unres/COMMON.SBRIDGE b/source/wham/src-M/include_unres/COMMON.SBRIDGE index 144a864..1380548 100644 --- a/source/wham/src-M/include_unres/COMMON.SBRIDGE +++ b/source/wham/src-M/include_unres/COMMON.SBRIDGE @@ -12,7 +12,7 @@ common /links_split/ link_start,link_end double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss logical dyn_ss,dyn_ss_mask - common /dyn_ssbond/ dtriss,atriss,btriss,ctriss, + common /dyn_ssbond/ dtriss,atriss,btriss,ctriss,Ht, & dyn_ssbond_ij(maxres,maxres), & idssb(maxdim),jdssb(maxdim), - & Ht,dyn_ss,dyn_ss_mask(maxres) + & dyn_ss,dyn_ss_mask(maxres) diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index d50cd03..750f03e 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -35,6 +35,7 @@ C character*16 key integer iparm double precision ip,mp +C write (iout,*) "KURWA" C C Body C @@ -55,6 +56,14 @@ C Assign virtual-bond length write (iout,*) "iparm",iparm," myparm",myparm c If reading not own parameters, skip assignment + call reada(weightcard,"D0CM",d0cm,3.78d0) + call reada(weightcard,"AKCM",akcm,15.1d0) + call reada(weightcard,"AKTH",akth,11.0d0) + call reada(weightcard,"AKCT",akct,12.0d0) + call reada(weightcard,"V1SS",v1ss,-1.08d0) + call reada(weightcard,"V2SS",v2ss,7.61d0) + call reada(weightcard,"V3SS",v3ss,13.7d0) + call reada(weightcard,"EBR",ebr,-5.50D0) call reada(controlcard,"DTRISS",dtriss,1.0D0) call reada(controlcard,"ATRISS",atriss,0.3D0) call reada(controlcard,"BTRISS",btriss,0.02D0) @@ -63,24 +72,45 @@ c If reading not own parameters, skip assignment C do i=1,maxres C dyn_ss_mask(i)=.false. C enddo +C ebr=-12.0D0 +c +c Old arbitrary potential - commented out. +c +c dbr= 4.20D0 +c fbr= 3.30D0 +c +c Constants of the disulfide-bond potential determined based on the RHF/6-31G** +c energy surface of diethyl disulfide. +c A. Liwo and U. Kozlowska, 11/24/03 +c + D0CM = 3.78d0 + AKCM = 15.1d0 + AKTH = 11.0d0 + AKCT = 12.0d0 + V1SS =-1.08d0 + V2SS = 7.61d0 + V3SS = 13.7d0 + do i=1,maxres-1 do j=i+1,maxres dyn_ssbond_ij(i,j)=1.0d300 enddo enddo call reada(controlcard,"HT",Ht,0.0D0) - if (dyn_ss) then - ss_depth=ebr/wsc-0.25*eps(1,1) - Ht=Ht/wsc-0.25*eps(1,1) - akcm=akcm*wstrain/wsc - akth=akth*wstrain/wsc - akct=akct*wstrain/wsc - v1ss=v1ss*wstrain/wsc - v2ss=v2ss*wstrain/wsc - v3ss=v3ss*wstrain/wsc - else - ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain - endif +C if (dyn_ss) then +C ss_depth=ebr/wsc-0.25*eps(1,1) +C write(iout,*) HT,wsc,eps(1,1),'KURWA' +C Ht=Ht/wsc-0.25*eps(1,1) + +C akcm=akcm*whpb/wsc +C akth=akth*whpb/wsc +C akct=akct*whpb/wsc +C v1ss=v1ss*whpb/wsc +C v2ss=v2ss*whpb/wsc +C v3ss=v3ss*whpb/wsc +C else +C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb +C endif if (iparm.eq.myparm .or. .not.separate_parset) then @@ -104,7 +134,7 @@ c wvdwpp=ww(16) wbond=ww(18) wsccor=ww(19) - + whpb=ww(15) endif call card_concat(controlcard,.false.) @@ -1217,7 +1247,7 @@ C C C Define the constants of the disulfide bridge C - ebr=-5.50D0 +C ebr=-12.0D0 c c Old arbitrary potential - commented out. c @@ -1228,21 +1258,36 @@ c Constants of the disulfide-bond potential determined based on the RHF/6-31G** c energy surface of diethyl disulfide. c A. Liwo and U. Kozlowska, 11/24/03 c - D0CM = 3.78d0 - AKCM = 15.1d0 - AKTH = 11.0d0 - AKCT = 12.0d0 - V1SS =-1.08d0 - V2SS = 7.61d0 - V3SS = 13.7d0 +C D0CM = 3.78d0 +C AKCM = 15.1d0 +C AKTH = 11.0d0 +C AKCT = 12.0d0 +C V1SS =-1.08d0 +C V2SS = 7.61d0 +C V3SS = 13.7d0 - if (lprint) then + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) +C write(iout,*) akcm,whpb,wsc,'KURWA' + Ht=Ht/wsc-0.25*eps(1,1) + + akcm=akcm*whpb/wsc + akth=akth*whpb/wsc + akct=akct*whpb/wsc + v1ss=v1ss*whpb/wsc + v2ss=v2ss*whpb/wsc + v3ss=v3ss*whpb/wsc + else + ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb + endif + +C if (lprint) then write (iout,'(/a)') "Disulfide bridge parameters:" write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, & ' v3ss:',v3ss - endif +C endif return end diff --git a/source/wham/src-M/ssMD.F b/source/wham/src-M/ssMD.F index 283adf3..5080b18 100644 --- a/source/wham/src-M/ssMD.F +++ b/source/wham/src-M/ssMD.F @@ -251,6 +251,7 @@ c-------END TESTING CODE e1=fac*fac*aa(itypi,itypj) e2=fac*bb(itypi,itypj) eij=eps1*eps2rt*eps3rt*(e1+e2) +C write(iout,*) eij,'TU?1' eps2der=eij*eps3rt eps3der=eij*eps2rt eij=eij*eps2rt*eps3rt @@ -267,7 +268,7 @@ c-------END TESTING CODE havebond=.true. ssd=rij-ssXs eij=ssA*ssd*ssd+ssB*ssd+ssC - +C write(iout,*) 'TU?2',ssc,ssd ed=2*akcm*ssd+akct*deltat12 pom1=akct*ssd pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi @@ -303,6 +304,7 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE h1=h_base(f1,hd1) h2=h_base(f2,hd2) eij=ssm*h1+Ht*h2 +C write(iout,*) eij,'TU?3' delta_inv=1.0d0/(xm-ssxm) deltasq_inv=delta_inv*delta_inv fac=ssm*hd1-Ht*hd2 @@ -325,6 +327,7 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE h1=h_base(f1,hd1) h2=h_base(f2,hd2) eij=Ht*h1+ljm*h2 +C write(iout,*) 'TU?4',ssA delta_inv=1.0d0/(ljxm-xm) deltasq_inv=delta_inv*delta_inv fac=Ht*hd1-ljm*hd2 @@ -396,7 +399,7 @@ c$$$ if (ed.gt.0.0d0) havebond=.true. c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE endif - + write(iout,*) 'havebond',havebond if (havebond) then #ifndef CLUST #ifndef WHAM -- 1.7.9.5