From: Adam Sieradzan Date: Thu, 15 Oct 2015 11:45:29 +0000 (+0200) Subject: small fix for Gly as 1st residue and in dynamic disulfides X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=d49e46757e9d141bbeff7e406cbc839fe6ab3952 small fix for Gly as 1st residue and in dynamic disulfides --- diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index cba6b5e..be6ed2a 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -23,7 +23,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' double precision fact(6) -cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot +Cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot cd print *,'nnt=',nnt,' nct=',nct C C Compute the side-chain and electrostatic interaction energy @@ -825,6 +825,7 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e +C write(iout,*) i,"i",iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -891,6 +892,8 @@ C & 'evdw',i,j,evdwij,'tss',evdw,evdw_t enddo! k ELSE ind=ind+1 +C write(iout,*) j,"j",istart(i,iint),iend(i,iint) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) @@ -1955,7 +1958,7 @@ C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ -cd write(iout,*) 'In EELEC' + write(iout,*) 'In EELEC' cd do i=1,nloctyp cd write(iout,*) 'Type',i cd write(iout,*) 'B1',B1(:,i) @@ -2013,15 +2016,16 @@ C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e gcorr_loc(i)=0.0d0 enddo do i=iatel_s,iatel_e - if (i.eq.1) then - if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1) cycle - else + write (iout,*) i,"i2",itype(i) + if (i.eq.1) cycle +C if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 +C & .or. itype(i+2).eq.ntyp1) cycle +C else if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 & .or. itype(i-1).eq.ntyp1 &) cycle - endif +C endif if (itel(i).eq.0) goto 1215 dxi=dc(1,i) dyi=dc(2,i) @@ -2041,6 +2045,8 @@ C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e num_conti=0 C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) +C write(iout,*) j,"j2" + if (j.eq.1) then if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 & .or.itype(j+2).eq.ntyp1 @@ -2051,6 +2057,7 @@ C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) & .or.itype(j-1).eq.ntyp1 &) cycle endif +C write(iout,*) j,"j2" C C) cycle if (itel(j).eq.0) goto 1216 @@ -2844,7 +2851,7 @@ C Cartesian derivatives & +0.5d0*(pizda(1,1)+pizda(2,2)) enddo endif - else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then + else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1.and.(i.gt.1)) then if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds & .or.((i+5).gt.nres) @@ -3682,7 +3689,7 @@ c & rad2deg*phii,rad2deg*phii1,ethetai c 1215 continue enddo ethetacnstr=0.0d0 -C print *,ithetaconstr_start,ithetaconstr_end,"TU" + print *,ntheta_constr,"TU" do i=1,ntheta_constr itheta=itheta_constr(i) thetiii=theta(itheta) @@ -3706,7 +3713,7 @@ C & i,itheta,rad2deg*thetiii, C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, C & gloc(itheta+nphi-2,icg) - endif +C endif enddo C Ufff.... We've done all this!!! return diff --git a/source/wham/src-M/ssMD.F b/source/wham/src-M/ssMD.F index 5080b18..bbf3206 100644 --- a/source/wham/src-M/ssMD.F +++ b/source/wham/src-M/ssMD.F @@ -150,11 +150,113 @@ c-------END TESTING CODE dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) dsci_inv=vbld_inv(i+nres) + xi=c(1,nres+i) + 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 + xj=c(1,nres+j) + 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 itypj=itype(j) - xj=c(1,nres+j)-c(1,nres+i) - yj=c(2,nres+j)-c(2,nres+i) - zj=c(3,nres+j)-c(3,nres+i) +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 + +C xj=c(1,nres+j)-c(1,nres+i) +C yj=c(2,nres+j)-c(2,nres+i) +C zj=c(3,nres+j)-c(3,nres+i) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -187,9 +289,9 @@ c om12=dxi*dxj+dyi*dyj+dzi*dzj ljXs=sig-sig0ij ljA=eps1*eps2rt**2*eps3rt**2 - ljB=ljA*bb(itypi,itypj) - ljA=ljA*aa(itypi,itypj) - ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) + ljB=ljA*bb + ljA=ljA*aa + ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0) ssXs=d0cm deltat1=1.0d0-om1 @@ -223,7 +325,7 @@ c-------TESTING CODE c Stop and plot energy and derivative as a function of distance if (checkstop) then ssm=ssC-0.25D0*ssB*ssB/ssA - ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) + ljm=-0.25D0*ljB*bb/aa if (ssm.lt.ljm .and. & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then nicheck=1000 @@ -248,8 +350,8 @@ c-------END TESTING CODE havebond=.false. ljd=rij-ljXs fac=(1.0D0/ljd)**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb eij=eps1*eps2rt*eps3rt*(e1+e2) C write(iout,*) eij,'TU?1' eps2der=eij*eps3rt @@ -316,8 +418,8 @@ C write(iout,*) eij,'TU?3' eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3) else havebond=.false. - ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) - d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB + ljm=-0.25D0*ljB*bb/aa + d_ljm(1)=-0.5D0*bb/aa*ljB d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt) d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt- + alf12/eps3rt)