X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-HCD%2Fenergy_p_new.F;h=db2e0438672a83312cffcdfd641d17da29c3de8a;hb=3d04f1babd4ae90f9a0b8f256a2ce0a2c5751f33;hp=5cc851c69766c8c8ab6b4e984967e004da84d3cf;hpb=1346fb3959c2eb0a370b11bc6ccad5e4cca27ec9;p=unres.git diff --git a/source/cluster/wham/src-HCD/energy_p_new.F b/source/cluster/wham/src-HCD/energy_p_new.F index 5cc851c..db2e043 100644 --- a/source/cluster/wham/src-HCD/energy_p_new.F +++ b/source/cluster/wham/src-HCD/energy_p_new.F @@ -180,14 +180,23 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t c write(iout,*) "TEST_ENE1 ehomology_constr=",ehomology_constr #ifdef DFA C BARTEK for dfa test! + edfadis=0.0d0 if (wdfa_dist.gt.0) call edfad(edfadis) c write(iout,*)'edfad is finished!', wdfa_dist,edfadis + edfator=0.0d0 if (wdfa_tor.gt.0) call edfat(edfator) c write(iout,*)'edfat is finished!', wdfa_tor,edfator + edfanei=0.0d0 if (wdfa_nei.gt.0) call edfan(edfanei) c write(iout,*)'edfan is finished!', wdfa_nei,edfanei + edfabet=0.0d0 if (wdfa_beta.gt.0) call edfab(edfabet) c write(iout,*)'edfab is finished!', wdfa_beta,edfabet +#else + edfadis=0.0d0 + edfator=0.0d0 + edfanei=0.0d0 + edfabet=0.0d0 #endif #ifdef SPLITELE @@ -511,6 +520,9 @@ C Bartek edfator = energia(29) edfanei = energia(30) edfabet = energia(31) + Eafmforc=0.0d0 + etube=0.0d0 + Uconst=0.0d0 #ifdef SPLITELE write(iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,wvdwpp, & estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1), @@ -674,6 +686,7 @@ cROZNICA xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) C Change 12/1/95 num_conti=0 C @@ -688,6 +701,10 @@ cd & 'iend=',iend(i,iint) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) C Change 12/1/95 to calculate four-body interactions rij=xj*xj+yj*yj+zj*zj rrij=1.0D0/rij @@ -857,6 +874,7 @@ c enddo xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) C C Calculate SC interaction energy. C @@ -867,6 +885,10 @@ C xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm @@ -972,6 +994,7 @@ c endif xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1004,9 +1027,13 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 - xj=c(1,nres+j)-xi - yj=c(2,nres+j)-yi - zj=c(3,nres+j)-zi + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1110,35 +1137,8 @@ c if (icall.gt.0) lprn=.true. yi=c(2,nres+i) zi=c(3,nres+i) C returning the ith atom to box - xi=mod(xi,boxxsize) - if (xi.lt.0) xi=xi+boxxsize - yi=mod(yi,boxysize) - if (yi.lt.0) yi=yi+boxysize - zi=mod(zi,boxzsize) - if (zi.lt.0) zi=zi+boxzsize - if ((zi.gt.bordlipbot) - &.and.(zi.lt.bordliptop)) then -C the energy transfer exist - if (zi.lt.buflipbot) then -C what fraction I am in - fracinbuf=1.0d0- - & ((zi-bordlipbot)/lipbufthick) -C lipbufthick is thickenes of lipid buffore - sslipi=sscalelip(fracinbuf) - ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) - sslipi=sscalelip(fracinbuf) - ssgradlipi=sscagradlip(fracinbuf)/lipbufthick - else - sslipi=1.0d0 - ssgradlipi=0.0 - endif - else - sslipi=0.0d0 - ssgradlipi=0.0 - endif - + call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1196,80 +1196,15 @@ c alf12=0.0D0 yj=c(2,nres+j) zj=c(3,nres+j) C returning jth atom to box - 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 - if ((zj.gt.bordlipbot) - &.and.(zj.lt.bordliptop)) then -C the energy transfer exist - if (zj.lt.buflipbot) then -C what fraction I am in - fracinbuf=1.0d0- - & ((zj-bordlipbot)/lipbufthick) -C lipbufthick is thickenes of lipid buffore - sslipj=sscalelip(fracinbuf) - ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zj.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) - sslipj=sscalelip(fracinbuf) - ssgradlipj=sscagradlip(fracinbuf)/lipbufthick - else - sslipj=1.0d0 - ssgradlipj=0.0 - endif - else - sslipj=0.0d0 - ssgradlipj=0.0 - endif - 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 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 -C if (aa.ne.aa_aq(itypi,itypj)) then - -C write(iout,*) "tu,", i,j,aa_aq(itypi,itypj)-aa, -C & bb_aq(itypi,itypj)-bb, -C & sslipi,sslipj -C endif - -C write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj) -C checking the distance - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 -C finding the closest - 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-xi)**2+(yj-yi)**2+(zj-zi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - subchap=1 - endif - enddo - enddo - enddo - if (subchap.eq.1) then - xj=xj_temp-xi - yj=yj_temp-yi - zj=zj_temp-zi - else - xj=xj_safe-xi - yj=yj_safe-yi - zj=zj_safe-zi - endif - + call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1391,6 +1326,8 @@ c if (icall.gt.0) lprn=.true. xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1425,9 +1362,15 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 - xj=c(1,nres+j)-xi - yj=c(2,nres+j)-yi - zj=c(3,nres+j)-zi + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) + call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1775,13 +1718,17 @@ C to calculate the el-loc multibody terms of various order. C c write(iout,*) 'SET_MATRICES nphi=',nphi,nres do i=3,nres+1 - if (i.gt. nnt+2 .and. i.lt.nct+2) then + ii=ireschain(i-2) + if (ii.eq.0) cycle + innt=chain_border(1,ii) + inct=chain_border(2,ii) + if (i.gt. innt+2 .and. i.lt.inct+2) then iti = itype2loc(itype(i-2)) else iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then - if (i.gt. nnt+1 .and. i.lt.nct+1) then + if (i.gt. innt+1 .and. i.lt.inct+1) then iti1 = itype2loc(itype(i-1)) else iti1=nloctyp @@ -1856,6 +1803,19 @@ c b2tilde(2,i-2)=-b2(2,i-2) write (iout,*) 'theta=', theta(i-1) #endif #else + if (i.gt. innt+2 .and. i.lt.inct+2) then +c if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itype2loc(itype(i-2)) + else + iti=nloctyp + endif +c write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti +c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itype2loc(itype(i-1)) + else + iti1=nloctyp + endif c if (i.gt. nnt+2 .and. i.lt.nct+2) then c iti = itype2loc(itype(i-2)) c else @@ -2239,12 +2199,7 @@ c end if xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*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) num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -2274,38 +2229,7 @@ c & .or. itype(i-1).eq.ntyp1 xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi -C Return atom into box, boxxsize is size of box in x dimension -c 194 continue -c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize -c if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize -C Condition for being inside the proper box -c if ((xmedi.gt.((0.5d0)*boxxsize)).or. -c & (xmedi.lt.((-0.5d0)*boxxsize))) then -c go to 194 -c endif -c 195 continue -c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize -c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize -C Condition for being inside the proper box -c if ((ymedi.gt.((0.5d0)*boxysize)).or. -c & (ymedi.lt.((-0.5d0)*boxysize))) then -c go to 195 -c endif -c 196 continue -c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize -c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize -C Condition for being inside the proper box -c if ((zmedi.gt.((0.5d0)*boxzsize)).or. -c & (zmedi.lt.((-0.5d0)*boxzsize))) then -c go to 196 -c endif - 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) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -2345,43 +2269,7 @@ c & .or. itype(i-1).eq.ntyp1 xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*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 -C xmedi=xmedi+xshift*boxxsize -C ymedi=ymedi+yshift*boxysize -C zmedi=zmedi+zshift*boxzsize - -C Return tom into box, boxxsize is size of box in x dimension -c 164 continue -c if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize -c if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize -C Condition for being inside the proper box -c if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or. -c & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then -c go to 164 -c endif -c 165 continue -c if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize -c if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize -C Condition for being inside the proper box -c if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or. -c & (ymedi.lt.((yshift-0.5d0)*boxysize))) then -c go to 165 -c endif -c 166 continue -c if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize -c if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize -cC Condition for being inside the proper box -c if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or. -c & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then -c go to 166 -c endif - -c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) + call to_box(xmedi,ymedi,zmedi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -2486,75 +2374,11 @@ C zj=c(3,j)+0.5D0*dzj-zmedi xj=c(1,j)+0.5D0*dxj yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*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 - if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ" - 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 if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC -c 174 continue -c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize -c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize -C Condition for being inside the proper box -c if ((xj.gt.((0.5d0)*boxxsize)).or. -c & (xj.lt.((-0.5d0)*boxxsize))) then -c go to 174 -c endif -c 175 continue -c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize -c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize -C Condition for being inside the proper box -c if ((yj.gt.((0.5d0)*boxysize)).or. -c & (yj.lt.((-0.5d0)*boxysize))) then -c go to 175 -c endif -c 176 continue -c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize -c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize -C Condition for being inside the proper box -c if ((zj.gt.((0.5d0)*boxzsize)).or. -c & (zj.lt.((-0.5d0)*boxzsize))) then -c go to 176 -c endif -C endif !endPBC condintion -C xj=xj-xmedi -C yj=yj-ymedi -C zj=zj-zmedi + call to_box(xj,yj,zj) + xj=boxshift(xj-xmedi,boxxsize) + yj=boxshift(yj-ymedi,boxysize) + zj=boxshift(zj-zmedi,boxzsize) rij=xj*xj+yj*yj+zj*zj - sss=sscale(sqrt(rij)) if (sss.eq.0.0d0) return sssgrad=sscagrad(sqrt(rij)) @@ -3075,7 +2899,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij c if (eel_loc_ij.ne.0) @@ -4038,12 +3862,7 @@ c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) C Returning the ith atom to box - xi=mod(xi,boxxsize) - if (xi.lt.0) xi=xi+boxxsize - yi=mod(yi,boxysize) - if (yi.lt.0) yi=yi+boxysize - zi=mod(zi,boxzsize) - if (zi.lt.0) zi=zi+boxzsize + call to_box(xi,yi,zi) do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -4058,44 +3877,10 @@ C Uncomment following three lines for Ca-p interactions yj=c(2,j) zj=c(3,j) C returning the jth atom to box - 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-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 -C Finding the closest jth atom - 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-xi)**2+(yj-yi)**2+(zj-zi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - subchap=1 - endif - enddo - enddo - enddo - if (subchap.eq.1) then - xj=xj_temp-xi - yj=yj_temp-yi - zj=zj_temp-zi - else - xj=xj_safe-xi - yj=yj_safe-yi - zj=zj_safe-zi - endif + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) C sss is scaling function for smoothing the cutoff gradient otherwise C the gradient would not be continuouse @@ -4520,6 +4305,10 @@ c estr1=0.0d0 c write (iout,*) "distchainmax",distchainmax do i=nnt+1,nct +#ifdef FIVEDIAG + if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) cycle + diff = vbld(i)-vbldp0 +#else if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) C do j=1,3 @@ -4537,6 +4326,9 @@ C write(iout,*) i,diff diff = vbld(i)-vbldp0 c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff endif +#endif + if (energy_dec) write (iout,'(a7,i5,4f7.3)') + & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff estr=estr+diff*diff do j=1,3 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) @@ -4555,8 +4347,9 @@ c nbi=nbondterm(iti) if (nbi.eq.1) then diff=vbld(i+nres)-vbldsc0(1,iti) -C write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, -C & AKSC(1,iti),AKSC(1,iti)*diff*diff + if (energy_dec) write (iout,*) + & i,iti,vbld(i+nres),vbldsc0(1,iti),diff, + & AKSC(1,iti),AKSC(1,iti)*diff*diff estr=estr+0.5d0*AKSC(1,iti)*diff*diff do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) @@ -10178,16 +9971,20 @@ c enddo c min_odl=minval(distancek) - do kk=1,constr_homology - if(l_homo(kk,ii)) then - min_odl=distancek(kk) - exit - endif - enddo - do kk=1,constr_homology - if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + if (nexl.gt.0) then + min_odl=0.0d0 + else + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) & min_odl=distancek(kk) - enddo + enddo + endif c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij