From: Adam Sieradzan Date: Mon, 20 Jul 2015 08:39:38 +0000 (+0200) Subject: debuging WHAM for constrains DEBUG ON X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=dd6b6b032024dfa86d8ed7622a72194c45b97cce debuging WHAM for constrains DEBUG ON --- diff --git a/source/wham/src-M/DIMENSIONS b/source/wham/src-M/DIMENSIONS index 154f3df..e9dac01 100644 --- a/source/wham/src-M/DIMENSIONS +++ b/source/wham/src-M/DIMENSIONS @@ -13,8 +13,8 @@ C Max. number of coarse-grain processors c parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres -c parameter (maxres=250) - parameter (maxres=100) + parameter (maxres=130) +C parameter (maxres=100) C Appr. max. number of interaction sites integer maxres2 parameter (maxres2=2*maxres) diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 6e19387..f184eed 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -2947,10 +2947,13 @@ C include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' + include 'COMMON.CONTROL' + include 'COMMON.IOUNITS' dimension ggg(3) ehpb=0.0D0 cd print *,'edis: nhpb=',nhpb,' fbr=',fbr cd print *,'link_start=',link_start,' link_end=',link_end + write(iout,*) link_end, "link_end" if (link_end.eq.0) return do i=link_start,link_end C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a @@ -2969,21 +2972,27 @@ C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. C & iabs(itype(jjj)).eq.1) then - + write(iout,*) constr_dist,"const" if (.not.dyn_ss .and. i.le.nss) then if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij - endif + endif !ii.gt.neres else if (ii.gt.nres .and. jj.gt.nres) then c Restraints from contact prediction dd=dist(ii,jj) if (constr_dist.eq.11) then +C ehpb=ehpb+fordepth(i)**4.0d0 +C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) ehpb=ehpb+fordepth(i)**4.0d0 & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) fac=fordepth(i)**4.0d0 & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + write(iout,*) ehpb,"atu?" +C ehpb,"tu?" +C fac=fordepth(i)**4.0d0 +C & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd else if (dhpb1(i).gt.0.0d0) then ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) @@ -3002,8 +3011,8 @@ C C Evaluate gradient. C fac=waga*rdis/dd - endif - endif + endif !end dhpb1(i).gt.0 + endif !end const_dist=11 do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@ -3015,7 +3024,26 @@ C ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) enddo - else + else !ii.gt.nres + write(iout,*) "before" + dd=dist(ii,jj) + write(iout,*) "after",dd + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) +C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd +C print *,ehpb,"tu?" + write(iout,*) ehpb,"btu?",dd,dhpb(i),dhpb1(i),fordepth(i) + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) diff --git a/source/wham/src-M/gnmr1.f b/source/wham/src-M/gnmr1.f index 905e746..8bfc43a 100644 --- a/source/wham/src-M/gnmr1.f +++ b/source/wham/src-M/gnmr1.f @@ -41,3 +41,33 @@ c------------------------------------------------------------------------------- return end c--------------------------------------------------------------------------------- + double precision function rlornmr1(y,ymin,ymax,sigma) + implicit none + double precision y,ymin,ymax,sigma + double precision wykl /4.0d0/ + if (y.lt.ymin) then + rlornmr1=(ymin-y)**wykl/((ymin-y)**wykl+sigma**wykl) + else if (y.gt.ymax) then + rlornmr1=(y-ymax)**wykl/((y-ymax)**wykl+sigma**wykl) + else + rlornmr1=0.0d0 + endif + return + end +c------------------------------------------------------------------------------ + double precision function rlornmr1prim(y,ymin,ymax,sigma) + implicit none + double precision y,ymin,ymax,sigma + double precision wykl /4.0d0/ + if (y.lt.ymin) then + rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/ + & ((ymin-y)**wykl+sigma**wykl)**2 + else if (y.gt.ymax) then + rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ + & ((y-ymax)**wykl+sigma**wykl)**2 + else + rlornmr1prim=0.0d0 + endif + return + end + diff --git a/source/wham/src-M/include_unres/COMMON.SBRIDGE b/source/wham/src-M/include_unres/COMMON.SBRIDGE index bea57e4..028f9ae 100644 --- a/source/wham/src-M/include_unres/COMMON.SBRIDGE +++ b/source/wham/src-M/include_unres/COMMON.SBRIDGE @@ -2,9 +2,10 @@ integer ns,nss,nfree,iss common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss, & ns,nss,nfree,iss(maxss) - double precision dhpb,dhpb1,forcon + double precision dhpb,dhpb1,forcon,fordepth integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim), + & fordepth(maxdim), & ihpb(maxdim),jhpb(maxdim),nhpb double precision weidis common /restraints/ weidis diff --git a/source/wham/src-M/read_dist_constr.F b/source/wham/src-M/read_dist_constr.F index 2e9ba3b..a97f797 100644 --- a/source/wham/src-M/read_dist_constr.F +++ b/source/wham/src-M/read_dist_constr.F @@ -13,9 +13,12 @@ character*500 controlcard logical lprn /.true./ write (iout,*) "Calling read_dist_constr" - write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup +C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup +C call flush(iout) + write(iout,*) "TU sie wywalam?" + call card_concat(controlcard,.false.) + write (iout,*) controlcard call flush(iout) - call card_concat(controlcard) call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NPAIR",npair_,0) call readi(controlcard,"NDIST",ndist_,0) @@ -100,16 +103,20 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) +C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", +C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) else read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + endif if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) endif - endif +C endif enddo + call hpb_partition call flush(iout) return end diff --git a/source/wham/src-M/wham_multparm.F b/source/wham/src-M/wham_multparm.F index e80a39d..de73250 100644 --- a/source/wham/src-M/wham_multparm.F +++ b/source/wham/src-M/wham_multparm.F @@ -58,6 +58,7 @@ c NaNQ initialization call flush(iout) call molread(*10) call flush(iout) + if (constr_dist.gt.0) call read_dist_constr #ifdef MPI write (iout,*) "Calling proc_groups" call proc_groups @@ -88,8 +89,8 @@ c NaNQ initialization call read_ref_structure(*10) call proc_cont call fragment_list - if (constr_dist.gt.0) call read_dist_constr endif +C if (constr_dist.gt.0) call read_dist_constr write (iout,*) "Begin read_database" call flush(iout) call read_database(*10)