X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fsc_move.F;h=c552ee06068b36564aa2c6568b1040f3fad88897;hb=325eda160c9ad2982501e091ca40606a29043712;hp=b6837fde271a4f35fe8a920d4432dcfdb3d3b80b;hpb=af72f8e89a5d33f0d86ba898d6c5bbbda4b25b84;p=unres.git diff --git a/source/unres/src_MD-M/sc_move.F b/source/unres/src_MD-M/sc_move.F index b6837fd..c552ee0 100644 --- a/source/unres/src_MD-M/sc_move.F +++ b/source/unres/src_MD-M/sc_move.F @@ -213,7 +213,7 @@ c Define what is meant by "neighbouring side-chain" c Don't do glycine or ends i=itype(res_pick) - if (i.eq.10 .or. i.eq.21) return + if (i.eq.10 .or. i.eq.ntyp1) return c Freeze everything (later will relax only selected side-chains) mask_r=.true. @@ -255,7 +255,7 @@ cd print *,'new ',(energy(k),k=0,n_ene) n_try=0 do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop) c Move the selected residue (don't worry if it fails) - call gen_side(itype(res_pick),theta(res_pick+1), + call gen_side(iabs(itype(res_pick)),theta(res_pick+1), + alph(res_pick),omeg(res_pick),fail) c Minimize the side-chains starting from the new arrangement @@ -719,11 +719,40 @@ c if (icall.eq.0) lprn=.true. do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + 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- + & ((positi-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-positi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -735,7 +764,7 @@ C do j=istart(i,iint),iend(i,iint) IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) dscj_inv=dsc_inv(itypj) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) @@ -757,14 +786,85 @@ 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 +C xj=c(1,nres+j)-xi +C yj=c(2,nres+j)-yi +C zj=c(3,nres+j)-zi + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) + 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- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/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 + + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=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-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 + dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) + sss=sscale((1.0d0/rij)/sigma(itypi,itypj)) + sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) + C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -783,16 +883,16 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') cd & restyp(itypi),i,restyp(itypj),j, cd & epsi,sigm,chi1,chi2,chip1,chip2, @@ -809,10 +909,13 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac + fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=ssgradlipi*evdwij + gg_lipj(3)=ssgradlipj*evdwij C Calculate angular part of the gradient. call sc_grad ENDIF