From: Cezary Czaplewski Date: Fri, 22 Jun 2018 12:01:10 +0000 (+0200) Subject: Adam's new constr_dist X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=8990ec14fa8e20cd2830e9e482631510fd058dcd Adam's new constr_dist --- diff --git a/source/unres/src_MD-M/COMMON.SBRIDGE b/source/unres/src_MD-M/COMMON.SBRIDGE index f866aa7..469add2 100644 --- a/source/unres/src_MD-M/COMMON.SBRIDGE +++ b/source/unres/src_MD-M/COMMON.SBRIDGE @@ -2,16 +2,18 @@ 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 - integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb + double precision dhpb,dhpb1,forcon,fordepth,xlscore,wboltzd + integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb,irestr_type common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim), - & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb + & fordepth(maxdim),xlscore(maxdim),wboltzd, + & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),irestr_type(maxdim), + & nhpb double precision weidis common /restraints/ weidis integer link_start,link_end common /links_split/ link_start,link_end - double precision Ht,dyn_ssbond_ij + double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss logical dyn_ss,dyn_ss_mask common /dyn_ssbond/ dyn_ssbond_ij(maxres,maxres), - & idssb(maxdim),jdssb(maxdim), - & Ht,dyn_ss,dyn_ss_mask(maxres) + & Ht,dtriss,atriss,btriss,ctriss,dyn_ss,dyn_ss_mask(maxres), + & idssb(maxdim),jdssb(maxdim) diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index ea2fe9e..8932003 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -195,8 +195,8 @@ C C Calculate the virtual-bond-angle energy. C if (wang.gt.0d0) then - call ebend(ebe) - else + call ebend(ebe) + else ebe=0 endif c print *,"Processor",myrank," computed UB" @@ -215,7 +215,7 @@ cd print *,'nterm=',nterm else etors=0 edihcnstr=0 - endif + endif if (constr_homology.ge.1) then call e_modeller(ehomology_constr) @@ -233,9 +233,9 @@ C C 6/23/01 Calculate double-torsional energy C if (wtor_d.gt.0) then - call etor_d(etors_d) + call etor_d(etors_d) else - etors_d=0 + etors_d=0 endif c print *,"Processor",myrank," computed Utord" C @@ -277,14 +277,14 @@ c write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr c write (iout,*) "From EsaxsC: Esaxs_constr",Esaxs_constr else Esaxs_constr = 0.0d0 - endif + endif c print *,"Processor",myrank," computed Ucorr" C C If performing constraint dynamics, call the constraint energy C after the equilibration time if(usampl.and.totT.gt.eq_time) then call EconstrQ - call Econstr_back + call Econstr_back else Uconst=0.0d0 Uconst_back=0.0d0 @@ -527,7 +527,7 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef DEBUG write (iout,*) "sum_gradient gsaxsc, gsaxsx" do i=0,nres - write (iout,'(i3,3e15.5,5x,3e15.5)') + write (iout,'(i3,3e15.5,5x,3e15.5)') & i,(gsaxsc(j,i),j=1,3),(gsaxsx(j,i),j=1,3) enddo call flush(iout) @@ -738,7 +738,7 @@ c enddo gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & 0.5d0*(wscp*gvdwc_scpp(j,i)+ - & welec*gelc_long(j,i) + + & welec*gelc_long(j,i)+ & wel_loc*gel_loc_long(j,i)+ & wcorr*gcorr_long(j,i)+ & wcorr5*gradcorr5_long(j,i)+ @@ -1530,6 +1530,7 @@ C include 'COMMON.SBRIDGE' logical lprn integer xshift,yshift,zshift + evdw=0.0D0 ccccc energy_dec=.false. C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -2779,7 +2780,7 @@ c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) c write(iout,*) 'b1=',b1(1,i-2) c write (iout,*) 'theta=', theta(i-1) - enddo + enddo #else if (i.gt. nnt+2 .and. i.lt.nct+2) then iti = itortyp(itype(i-2)) @@ -2796,10 +2797,10 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then b1(2,i-2)=b(5,iti) b2(1,i-2)=b(2,iti) b2(2,i-2)=b(4,iti) - b1tilde(1,i-2)=b1(1,i-2) - b1tilde(2,i-2)=-b1(2,i-2) - b2tilde(1,i-2)=b2(1,i-2) - b2tilde(2,i-2)=-b2(2,i-2) + b1tilde(1,i-2)= b1(1,i-2) + b1tilde(2,i-2)=-b1(2,i-2) + b2tilde(1,i-2)= b2(1,i-2) + b2tilde(2,i-2)=-b2(2,i-2) EE(1,2,i-2)=eeold(1,2,iti) EE(2,1,i-2)=eeold(2,1,iti) EE(2,2,i-2)=eeold(2,2,iti) @@ -2893,7 +2894,7 @@ cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) cd write (iout,*) 'b2',b2(:,iti) cd write (iout,*) "phi(",i,")=",phi(i)," sin1",sin1," cos1",cos1 -cd write (iout,*) 'Ug',Ug(:,:,i-2) +cd write (iout,*) 'Ug',Ug(:,:,i-2) c if (i .gt. iatel_s+2) then if (i .gt. nnt+2) then call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2)) @@ -3730,6 +3731,7 @@ C erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij + * * Radial derivatives. First process both termini of the fragment (i,j) * @@ -5189,10 +5191,16 @@ C include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 + do i=1,3 + ggg(i)=0.0d0 + enddo +C write (iout,*) ,"link_end",link_end,constr_dist cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr -cd write(iout,*)'link_start=',link_start,' link_end=',link_end +c write(iout,*)'link_start=',link_start,' link_end=',link_end, +c & " constr_dist",constr_dist 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 @@ -5217,38 +5225,75 @@ cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds if (.not.dyn_ss .and. i.le.nss) then C 15/02/13 CC dynamic SSbond - additional check - if (ii.gt.nres - & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then - call ssbond_ene(iii,jjj,eij) - ehpb=ehpb+2*eij + 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 cd write (iout,*) "eij",eij - else +cd & ' waga=',waga,' fac=',fac +! else if (ii.gt.nres .and. jj.gt.nres) then + else C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) + if (irestr_type(i).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 + if (energy_dec) write (iout,'(a6,2i5,6f10.3,i5)') + & "edisL",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i), + & ehpb,irestr_type(i) + else if (irestr_type(i).eq.10) then +c AL 6//19/2018 cross-link restraints + xdis = 0.5d0*(dd/forcon(i))**2 + expdis = dexp(-xdis) +c aux=(dhpb(i)+dhpb1(i)*xdis)*expdis+fordepth(i) + aux=(dhpb(i)+dhpb1(i)*xdis*xdis)*expdis+fordepth(i) +c write (iout,*)"HERE: xdis",xdis," expdis",expdis," aux",aux, +c & " wboltzd",wboltzd + ehpb=ehpb-wboltzd*xlscore(i)*dlog(aux) +c fac=-wboltzd*(dhpb1(i)*(1.0d0-xdis)-dhpb(i)) + fac=-wboltzd*xlscore(i)*(dhpb1(i)*(2.0d0-xdis)*xdis-dhpb(i)) + & *expdis/(aux*forcon(i)**2) + if (energy_dec) write(iout,'(a6,2i5,6f10.3,i5)') + & "edisX",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i), + & -wboltzd*xlscore(i)*dlog(aux),irestr_type(i) + else if (irestr_type(i).eq.2) then +c Quartic restraints + ehpb=ehpb+forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)') + & "edisQ",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i), + & forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)),irestr_type(i) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd + else +c Quadratic restraints rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+0.5d0*waga*rdis*rdis + if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)') + & "edisS",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i), + & 0.5d0*waga*rdis*rdis,irestr_type(i) C C Evaluate gradient. C fac=waga*rdis/dd -cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, -cd & ' waga=',waga,' fac=',fac - do j=1,3 - ggg(j)=fac*(c(j,jj)-c(j,ii)) - enddo + endif +c Calculate Cartesian gradient + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). if (iii.lt.ii) then - do j=1,3 - ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) - ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) - enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo endif cgrad do j=iii,jjj-1 cgrad do k=1,3 @@ -5261,7 +5306,6 @@ cgrad enddo enddo endif enddo - ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -5400,6 +5444,7 @@ C NO vbldp0 is the equlibrium lenght of spring for peptide group c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) c endif enddo + estr=0.5d0*AKP*estr+estr1 c c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included @@ -5580,6 +5625,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2. if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo + C Ufff.... We've done all this!!! return end @@ -5896,6 +5942,7 @@ c lprn1=.false. if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg) enddo + return end #endif @@ -8587,7 +8634,7 @@ c & ' eij',eij,' eesij',ees0pij,ees0mij,' and ',k,l c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' energy=',ekont*ees, c & 'gradcorr_long' C Calculate the multi-body contribution to energy. -c ecorr=ecorr+ekont*ees +C ecorr=ecorr+ekont*ees C Calculate multi-body contributions to the gradient. coeffpees0pij=coeffp*ees0pij coeffmees0mij=coeffm*ees0mij @@ -11082,11 +11129,11 @@ c write (2,*) "rr",rr else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then gamm=(rr-(r_cut-rlamb))/rlamb sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0) - else + else sscale2=0d0 endif - return - end + return + end C----------------------------------------------------------------------- double precision function sscalgrad2(r,r_cut,r0,rlamb) implicit none @@ -11101,11 +11148,11 @@ C----------------------------------------------------------------------- else sscalgrad2=-gamm*(6*gamm-6.0d0)/rlamb endif - else + else sscalgrad2=0.0d0 endif - return - end + return + end c---------------------------------------------------------------------------- subroutine e_saxs(Esaxs_constr) implicit none diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 7b3f9a8..c511fb5 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -41,12 +41,18 @@ C Print restraint information #ifdef MPI if (.not. out1file .or. me.eq.king) then #endif - if (nhpb.gt.nss) - &write (iout,'(a,i5,a)') "The following",nhpb-nss, - & " distance constraints have been imposed" - do i=nss+1,nhpb - write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i) - enddo + if (nhpb.gt.nss) then + write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)') + & "The following",nhpb-nss, + & " distance restraints have been imposed:", + & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V", + & " score"," type" + do i=nss+1,nhpb + write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i), + & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i), + & irestr_type(i) + enddo + endif #ifdef MPI endif #endif @@ -2545,18 +2551,40 @@ c------------------------------------------------------------------------------- integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard -c write (iout,*) "Calling read_dist_constr" + logical normalize,next + integer restr_type + double precision xlink(4,0:4) / +c a b c sigma + & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential + & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL + & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH + & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH + & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS +c print *, "WCHODZE" + write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) + next=.true. + + DO WHILE (next) + call card_concat(controlcard) + next = index(controlcard,"NEXT").gt.0 + call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist) + write (iout,*) "restr_type",restr_type + call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NPAIR",npair_,0) call readi(controlcard,"NDIST",ndist_,0) call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) + if (restr_type.eq.10) + & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0) call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0) call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0) call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0) call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0) + normalize = index(controlcard,"NORMALIZE").gt.0 + write (iout,*) "WBOLTZD",wboltzd c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ c write (iout,*) "IFRAG" c do i=1,nfrag_ @@ -2566,46 +2594,38 @@ c write (iout,*) "IPAIR" c do i=1,npair_ c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) c enddo - if (.not.refstr .and. nfrag.gt.0) then - write (iout,*) - & "ERROR: no reference structure to compute distance restraints" - write (iout,*) - & "Restraints must be specified explicitly (NDIST=number)" - stop - endif - if (nfrag.lt.2 .and. npair.gt.0) then - write (iout,*) "ERROR: Less than 2 fragments specified", - & " but distance restraints between pairs requested" - stop - endif - call flush(iout) + if (nfrag_.gt.0) write (iout,*) + & "Distance restraints as generated from reference structure" do i=1,nfrag_ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup if (ifrag_(2,i).gt.nstart_sup+nsup-1) & ifrag_(2,i)=nstart_sup+nsup-1 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) - call flush(iout) - if (wfrag_(i).gt.0.0d0) then +c call flush(iout) + if (wfrag_(i).eq.0.0d0) cycle do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) c write (iout,*) "j",j," k",k ddjk=dist(j,k) - if (constr_dist.eq.1) then - nhpb=nhpb+1 - ihpb(nhpb)=j - jhpb(nhpb)=k + if (restr_type.eq.1) then + nhpb=nhpb+1 + irestr_type(nhpb)=1 + ihpb(nhpb)=j + jhpb(nhpb)=k dhpb(nhpb)=ddjk - forcon(nhpb)=wfrag_(i) + forcon(nhpb)=wfrag_(i) else if (constr_dist.eq.2) then if (ddjk.le.dist_cut) then nhpb=nhpb+1 + irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) endif - else + else if (restr_type.eq.3) then nhpb=nhpb+1 + irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk @@ -2613,18 +2633,17 @@ c write (iout,*) "j",j," k",k endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", + & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", + write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo enddo - endif enddo do i=1,npair_ - if (wpair_(i).gt.0.0d0) then + if (wpair_(i).eq.0.0d0) cycle ii = ipair_(1,i) jj = ipair_(2,i) if (ii.gt.jj) then @@ -2634,39 +2653,161 @@ c write (iout,*) "j",j," k",k endif do j=ifrag_(1,ii),ifrag_(2,ii) do k=ifrag_(1,jj),ifrag_(2,jj) - nhpb=nhpb+1 - ihpb(nhpb)=j - jhpb(nhpb)=k - forcon(nhpb)=wpair_(i) - dhpb(nhpb)=dist(j,k) + if (restr_type.eq.1) then + nhpb=nhpb+1 + irestr_type(nhpb)=1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + else if (constr_dist.eq.2) then + if (ddjk.le.dist_cut) then + nhpb=nhpb+1 + irestr_type(nhpb)=1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + endif + else if (restr_type.eq.3) then + nhpb=nhpb+1 + irestr_type(nhpb)=1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) + endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo enddo - endif enddo + +c print *,ndist_ + write (iout,*) "Distance restraints as read from input" do i=1,ndist_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) - if (forcon(nhpb+1).gt.0.0d0) then + if (restr_type.eq.11) then + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), + & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) +c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) + if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle + nhpb=nhpb+1 + irestr_type(nhpb)=11 +#ifdef MPI + if (.not.out1file .or. me.eq.king) + & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), + & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) +#else + write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), + & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) +#endif + if (ibecarb(nhpb).gt.0) then + ihpb(nhpb)=ihpb(nhpb)+nres + jhpb(nhpb)=jhpb(nhpb)+nres + endif + else if (constr_dist.eq.10) then +c Cross-lonk Markov-like potential + call card_concat(controlcard) + call readi(controlcard,"ILINK",ihpb(nhpb+1),0) + call readi(controlcard,"JLINK",jhpb(nhpb+1),0) + ibecarb(nhpb+1)=0 + if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1 + if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle + if (index(controlcard,"ZL").gt.0) then + link_type=1 + else if (index(controlcard,"ADH").gt.0) then + link_type=2 + else if (index(controlcard,"PDH").gt.0) then + link_type=3 + else if (index(controlcard,"DSS").gt.0) then + link_type=4 + else + link_type=0 + endif + call reada(controlcard,"AXLINK",dhpb(nhpb+1), + & xlink(1,link_type)) + call reada(controlcard,"BXLINK",dhpb1(nhpb+1), + & xlink(2,link_type)) + call reada(controlcard,"CXLINK",fordepth(nhpb+1), + & xlink(3,link_type)) + call reada(controlcard,"SIGMA",forcon(nhpb+1), + & xlink(4,link_type)) + call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0) +c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1), +c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1) + if (forcon(nhpb+1).le.0.0d0 .or. + & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle + nhpb=nhpb+1 + irestr_type(nhpb)=10 + if (ibecarb(nhpb).gt.0) then + ihpb(nhpb)=ihpb(nhpb)+nres + jhpb(nhpb)=jhpb(nhpb)+nres + endif +#ifdef MPI + if (.not.out1file .or. me.eq.king) + & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), + & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), + & irestr_type(nhpb) +#else + write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), + & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), + & irestr_type(nhpb) +#endif + else +C print *,"in else" + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), + & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1) + if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 - dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (dhpb1(nhpb).eq.0.0d0) then + irestr_type(nhpb)=1 + else + irestr_type(nhpb)=2 + endif + if (ibecarb(nhpb).gt.0) then + ihpb(nhpb)=ihpb(nhpb)+nres + jhpb(nhpb)=jhpb(nhpb)+nres + endif + if (dhpb(nhpb).eq.0.0d0) + & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", - & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", - & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #endif endif +C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +C if (forcon(nhpb+1).gt.0.0d0) then +C nhpb=nhpb+1 +C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) enddo - call flush(iout) + + ENDDO ! next + + fordepthmax=0.0d0 + if (normalize) then + do i=nss+1,nhpb + if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax) + & fordepthmax=fordepth(i) + enddo + do i=nss+1,nhpb + if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax + enddo + endif return end c------------------------------------------------------------------------------- diff --git a/source/unres/src_MD-M/unres.F b/source/unres/src_MD-M/unres.F index 0e924f7..073c94f 100644 --- a/source/unres/src_MD-M/unres.F +++ b/source/unres/src_MD-M/unres.F @@ -303,7 +303,7 @@ crc overlap test if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) call briefout(0,etot) endif - if (outpdb) call pdbout(etot,titel(:32),ipdb) + if (outpdb) call pdbout(etot,titel(:50),ipdb) if (outmol2) call mol2out(etot,titel(:32)) return end @@ -341,7 +341,7 @@ c--------------------------------------------------------------------------- call enerprint(energy(0)) call intout call briefout(0,etot) - if (outpdb) call pdbout(etot,titel(:32),ipdb) + if (outpdb) call pdbout(etot,titel(:50),ipdb) if (outmol2) call mol2out(etot,titel(:32)) if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) write (iout,'(a,i3)') 'SUMSL return code:',iretcode