X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=cba6b5e3c1a1ca116c5c77ce4237a86b04b82a7f;hb=25618f9f83673a7063414fe1e17415d138f58da8;hp=fcc52c0cad5e74db4a9ce6558bb9f19720e7810c;hpb=59935311712ea011d55e8169709e7c8076cb0b75;p=unres.git diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index fcc52c0..cba6b5e 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -44,11 +44,13 @@ C Gay-Berne potential (shifted LJ, angular dependence). goto 106 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). 105 call egbv(evdw,evdw_t) +C write(iout,*) 'po elektostatyce' C C Calculate electrostatic (H-bonding) energy of the main chain. C - 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) -C + 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) +C write(iout,*) 'po eelec' + C Calculate excluded-volume interaction energy between peptide groups C and side chains. C @@ -56,8 +58,9 @@ C c c Calculate the bond-stretching energy c + call ebond(estr) -c write (iout,*) "estr",estr +C write (iout,*) "estr",estr C C Calculate the disulfide-bridge and other energy and the contributions C from other distance constraints. @@ -67,13 +70,14 @@ cd print *,'EHPB exitted succesfully.' C C Calculate the virtual-bond-angle energy. C - call ebend(ebe) +C print *,'Bend energy finished.' + call ebend(ebe,ethetacnstr) cd print *,'Bend energy finished.' C C Calculate the SC local energy. C call esc(escloc) -cd print *,'SCLOC energy finished.' +C print *,'SCLOC energy finished.' C C Calculate the virtual-bond torsional energy. C @@ -87,6 +91,11 @@ C C 21/5/07 Calculate local sicdechain correlation energy C call eback_sc_corr(esccor) + + if (wliptran.gt.0) then + call Eliptransfer(eliptran) + endif + C C 12/1/95 Multi-body terms C @@ -107,20 +116,22 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees & +wvdwpp*evdw1 & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr + & +wliptran*eliptran #else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr + & +wliptran*eliptran #endif energia(0)=etot energia(1)=evdw @@ -154,6 +165,8 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t + energia(24)=ethetacnstr + energia(22)=eliptran c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -192,10 +205,12 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(2)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) enddo #else do i=1,nct @@ -211,10 +226,12 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(1)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) enddo #endif enddo @@ -228,9 +245,11 @@ C & +wturn3*fact(2)*gel_loc_turn3(i) & +wturn6*fact(5)*gel_loc_turn6(i) & +wel_loc*fact(2)*gel_loc_loc(i) - & +wsccor*fact(1)*gsccor_loc(i) +c & +wsccor*fact(1)*gsccor_loc(i) +c BYLA ROZNICA Z CLUSTER< OSTATNIA LINIA DODANA enddo endif + if (dyn_ss) call dyn_set_nss return end C------------------------------------------------------------------------ @@ -268,6 +287,8 @@ C------------------------------------------------------------------------ esccor=energia(19) edihcnstr=energia(20) estr=energia(18) + ethetacnstr=energia(24) + eliptran=energia(22) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@ -276,7 +297,8 @@ C------------------------------------------------------------------------ & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), - & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot + & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss, + & eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -298,7 +320,9 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond, @@ -307,7 +331,7 @@ C------------------------------------------------------------------------ & ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2), & eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3), & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, - & edihcnstr,ebr*nss,etot + & edihcnstr,ethetacnstr,ebr*nss,eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -328,7 +352,9 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -360,17 +386,20 @@ C integer icant external icant cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon +c ROZNICA z cluster do i=1,210 do j=1,2 eneps_temp(j,i)=0.0d0 enddo enddo +cROZNICA + evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -383,8 +412,8 @@ C cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -394,19 +423,22 @@ C Change 12/1/95 to calculate four-body interactions c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=e1+e2 ij=icant(itypi,itypj) +c ROZNICA z cluster eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij) eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij +c + cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij @@ -541,9 +573,9 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -552,8 +584,8 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -564,8 +596,8 @@ C rij=1.0D0/r_inv_ij r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=e_augm+e1+e2 ij=icant(itypi,itypj) eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm) @@ -578,7 +610,7 @@ cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij @@ -654,9 +686,9 @@ c else c endif ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -670,8 +702,8 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) @@ -710,8 +742,8 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives fac=(rrij*sigsq)**expon2 - 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 @@ -721,15 +753,15 @@ C to its derivatives eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux & /dabs(eps(itypi,itypj)) eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij endif if (calc_grad) then 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 write (iout,'(2(a3,i3,2x),15(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, @@ -776,6 +808,7 @@ C include 'COMMON.ENEPS' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SBRIDGE' logical lprn common /srutu/icall integer icant @@ -792,12 +825,42 @@ 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 - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) 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 + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -807,9 +870,29 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij +C write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)') +C & 'evdw',i,j,evdwij,' ss',evdw,evdw_t +C triple bond artifac removal + do k=j+1,iend(i,iint) +C search over all next residues + if (dyn_ss_mask(k)) then +C check if they are cysteins +C write(iout,*) 'k=',k + call triple_ssbond_ene(i,j,k,evdwij) +C call the energy function that removes the artifical triple disulfide +C bond the soubroutine is located in ssMD.F + evdw=evdw+evdwij +C write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)') +C & 'evdw',i,j,evdwij,'tss',evdw,evdw_t + endif!dyn_ss_mask(k) + enddo! k + ELSE ind=ind+1 - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) @@ -831,17 +914,89 @@ 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) +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 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 + dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) c write (iout,*) i,j,xj,yj,zj 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)) + if (sss.le.0.0) cycle C Calculate angle-dependent terms of energy and contributions to their C derivatives. + call sc_angular sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) @@ -855,16 +1010,16 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! 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 - if (bb(itypi,itypj).gt.0) then - evdw=evdw+evdwij + if (bb.gt.0) then + evdw=evdw+evdwij*sss else - evdw_t=evdw_t+evdwij + evdw_t=evdw_t+evdwij*sss endif ij=icant(itypi,itypj) aux=eps1*eps2rt**2*eps3rt**2 @@ -875,8 +1030,8 @@ c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj, c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)), c & aux*e2/eps(itypi,itypj) c 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 #ifdef DEBUG write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, @@ -893,6 +1048,7 @@ 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 @@ -900,6 +1056,8 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad endif +C write(iout,*) "partial sum", evdw, evdw_t + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -942,9 +1100,9 @@ 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 - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -958,8 +1116,8 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) r0ij=r0(itypi,itypj) @@ -1005,15 +1163,15 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! 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 fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm evdwij=evdwij*eps2rt*eps3rt - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij+e_augm else evdw_t=evdw_t+evdwij+e_augm @@ -1724,6 +1882,8 @@ c print *,"itilde3 i iti iti1",i,iti,iti1 do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) enddo +C write (iout,*) 'mumu',i,b1(1,iti),Ub2(1,i-2) + C Vectors and matrices dependent on a single virtual-bond dihedral. call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1)) call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) @@ -1846,14 +2006,22 @@ cd enddo do i=1,nres num_cont_hb(i)=0 enddo -cd print '(a)','Enter EELEC' -cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e +C print '(a)','Enter EELEC' +C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e do i=1,nres gel_loc_loc(i)=0.0d0 gcorr_loc(i)=0.0d0 enddo do i=iatel_s,iatel_e - if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (i.eq.1) then + if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 + & .or. itype(i+2).eq.ntyp1) cycle + 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 if (itel(i).eq.0) goto 1215 dxi=dc(1,i) dyi=dc(2,i) @@ -1864,10 +2032,27 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e 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 num_conti=0 -c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) +C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) - if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle + if (j.eq.1) then + if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 + & .or.itype(j+2).eq.ntyp1 + &) cycle + else + if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 + & .or.itype(j+2).eq.ntyp1 + & .or.itype(j-1).eq.ntyp1 + &) cycle + endif +C +C) cycle if (itel(j).eq.0) goto 1216 ind=ind+1 iteli=itel(i) @@ -1889,10 +2074,49 @@ C End diagnostics dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - 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 + 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 rij=xj*xj+yj*yj+zj*zj + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) rrmij=1.0D0/rij rij=dsqrt(rij) rmij=1.0D0/rij @@ -1916,16 +2140,21 @@ c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) ees=ees+eesij - evdw1=evdw1+evdwij -cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') -cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, -cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, -cd & xmedi,ymedi,zmedi,xj,yj,zj + evdw1=evdw1+evdwij*sss +c write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') +c &'evdw1',i,j,evdwij +c &,iteli,itelj,aaa,evdw1 + +C write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij +c write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') +c & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, +c & 1.0D0/dsqrt(rrmij),evdwij,eesij, +c & xmedi,ymedi,zmedi,xj,yj,zj C C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE - facvdw=-6*rrmij*(ev1+evdwij) + facvdw=-6*rrmij*(ev1+evdwij)*sss facel=-3*rrmij*(el1+eesij) fac1=fac erij(1)=xj*rmij @@ -1951,9 +2180,18 @@ C gelc(l,k)=gelc(l,k)+ggg(l) enddo enddo - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj +C ggg(1)=facvdw*xj +C ggg(2)=facvdw*yj +C ggg(3)=facvdw*zj + if (sss.gt.0.0) then + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + else + ggg(1)=0.0 + ggg(2)=0.0 + ggg(3)=0.0 + endif do k=1,3 ghalf=0.5D0*ggg(k) gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -1968,7 +2206,7 @@ C enddo enddo #else - facvdw=ev1+evdwij + facvdw=(ev1+evdwij)*sss facel=el1+eesij fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) @@ -2266,8 +2504,11 @@ C Check the loc-el terms by numerical integration C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) -cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij -cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) +c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij +C write (iout,'(a6,2i5,0pf7.3)') +C & 'eelloc',i,j,eel_loc_ij +C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4) +c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) eel_loc=eel_loc+eel_loc_ij C Partial derivatives in virtual-bond dihedral angles gamma if (calc_grad) then @@ -2603,7 +2844,17 @@ C Cartesian derivatives & +0.5d0*(pizda(1,1)+pizda(2,2)) enddo endif - else if (j.eq.i+3 .and. itype(i+2).ne.21) then + else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) 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) + & .or.((i-1).le.0) +C end of changes suggested by Ana + & .or. itype(i+3).eq.ntyp1 + & .or. itype(i+4).eq.ntyp1 + & .or. itype(i+5).eq.ntyp1 + & .or. itype(i).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1) goto 178 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Fourth-order contributions @@ -2743,6 +2994,7 @@ C Remaining derivatives of this turn contribution gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) enddo endif + 178 continue endif return end @@ -2804,7 +3056,7 @@ cd print '(a)','Enter ESCP' c write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e, c & ' scal14',scal14 do i=iatscp_s,iatscp_e - if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) c write (iout,*) "i",i," iteli",iteli," nscp_gr",nscp_gr(i), c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) @@ -2812,37 +3064,90 @@ c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) xi=0.5D0*(c(1,i)+c(1,i+1)) 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 do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + xj=c(1,j) + 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 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 + sss=sscale(1.0d0/(dsqrt(rrij))) + if (sss.le.0.0d0) cycle + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 - evdw2_14=evdw2_14+e1+e2 + evdw2_14=evdw2_14+(e1+e2)*sss endif evdwij=e1+e2 -c write (iout,*) i,j,evdwij - evdw2=evdw2+evdwij +c write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') +c & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), +c & bad(itypj,iteli) + evdw2=evdw2+evdwij*sss if (calc_grad) then C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C - fac=-(evdwij+e1)*rrij + fac=-(evdwij+e1)*rrij*sss + fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -2907,10 +3212,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 +C 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 @@ -2927,24 +3235,98 @@ C iii and jjj point to the residues for which the distance is assigned. endif C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +C & iabs(itype(jjj)).eq.1) then +C 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 - else -C Calculate the distance between the two points and its difference from the -C target distance. - dd=dist(ii,jj) - rdis=dd-dhpb(i) + 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 +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd +C 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)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + 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 +c write (iout,*) "beta reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + 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 + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo + else !ii.gt.nres +C write(iout,*) "before" + dd=dist(ii,jj) +C 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?" +C write(iout,*) ehpb,"btu?", +C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i) +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd + 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) + waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis 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 + fac=waga*rdis/dd + endif + endif + do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@ -2964,7 +3346,7 @@ C Cartesian gradient in the SC vectors (ghpbx). enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -2987,7 +3369,7 @@ C include 'COMMON.VAR' include 'COMMON.IOUNITS' double precision erij(3),dcosom1(3),dcosom2(3),gg(3) - itypi=itype(i) + itypi=iabs(itype(i)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -2995,7 +3377,7 @@ C dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) dsci_inv=dsc_inv(itypi) - itypj=itype(j) + itypj=iabs(itype(j)) dscj_inv=dsc_inv(itypj) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -3073,39 +3455,45 @@ c logical energy_dec /.false./ double precision u(3),ud(3) estr=0.0d0 - write (iout,*) "distchainmax",distchainmax + estr1=0.0d0 +c write (iout,*) "distchainmax",distchainmax do i=nnt+1,nct - if (itype(i-1).eq.21 .or. itype(i).eq.21) then - estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) - do j=1,3 - gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) - & *dc(j,i-1)/vbld(i) - enddo - if (energy_dec) write(iout,*) - & "estr1",i,vbld(i),distchainmax, - & gnmr1(vbld(i),-1.0d0,distchainmax) - 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 +C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) +C & *dc(j,i-1)/vbld(i) +C enddo +C if (energy_dec) write(iout,*) +C & "estr1",i,vbld(i),distchainmax, +C & gnmr1(vbld(i),-1.0d0,distchainmax) +C else + if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then + diff = vbld(i)-vbldpDUM + else diff = vbld(i)-vbldp0 c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff + endif estr=estr+diff*diff do j=1,3 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) enddo - endif - +C endif +C write (iout,'(a7,i5,4f7.3)') +C & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff enddo - estr=0.5d0*AKP*estr + estr=0.5d0*AKP*estr+estr1 c c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included c do i=nnt,nct - iti=itype(i) - if (iti.ne.10 .and. iti.ne.21) then + iti=iabs(itype(i)) + if (iti.ne.10 .and. iti.ne.ntyp1) then 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 +C write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, +C & 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) @@ -3147,7 +3535,7 @@ c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -3164,27 +3552,48 @@ C include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.FFIELD' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it double precision y(2),z(2) delta=0.02d0*pi - time11=dexp(-2*time) - time12=1.0d0 +c time11=dexp(-2*time) +c time12=1.0d0 etheta=0.0D0 c write (iout,*) "nres",nres c write (*,'(a,i2)') 'EBEND ICG=',icg c write (iout,*) ithet_start,ithet_end do i=ithet_start,ithet_end - if (itype(i-1).eq.21) cycle +C if (itype(i-1).eq.ntyp1) cycle + if (i.le.2) cycle + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) it=itype(i-1) - if (i.gt.3 .and. itype(i-2).ne.21) then + ichir1=isign(1,itype(i-2)) + ichir2=isign(1,itype(i)) + if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1)) + if (itype(i).eq.10) ichir2=isign(1,itype(i-1)) + if (itype(i-1).eq.10) then + itype1=isign(10,itype(i-2)) + ichir11=isign(1,itype(i-2)) + ichir12=isign(1,itype(i-2)) + itype2=isign(10,itype(i)) + ichir21=isign(1,itype(i)) + ichir22=isign(1,itype(i)) + endif + if (i.eq.3) then + y(1)=0.0D0 + y(2)=0.0D0 + else + + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) - icrc=0 - call proc_proc(phii,icrc) +c icrc=0 +c call proc_proc(phii,icrc) if (icrc.eq.1) phii=150.0 #else phii=phi(i) @@ -3195,11 +3604,12 @@ C Zero the energy function and its derivative at 0 or pi. y(1)=0.0D0 y(2)=0.0D0 endif - if (i.lt.nres .and. itype(i).ne.21) then + endif + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) - icrc=0 - call proc_proc(phii1,icrc) +c icrc=0 +c call proc_proc(phii1,icrc) if (icrc.eq.1) phii1=150.0 phii1=pinorm(phii1) z(1)=cos(phii1) @@ -3217,8 +3627,12 @@ C dependent on the adjacent virtual-bond-valence angles (gamma1 & gamma2). C In following comments this theta will be referred to as t_c. thet_pred_mean=0.0d0 do k=1,2 - athetk=athet(k,it) - bthetk=bthet(k,it) + athetk=athet(k,it,ichir1,ichir2) + bthetk=bthet(k,it,ichir1,ichir2) + if (it.eq.10) then + athetk=athet(k,itype1,ichir11,ichir12) + bthetk=bthet(k,itype2,ichir21,ichir22) + endif thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) enddo c write (iout,*) "thet_pred_mean",thet_pred_mean @@ -3226,8 +3640,16 @@ c write (iout,*) "thet_pred_mean",thet_pred_mean thet_pred_mean=thet_pred_mean*ss+a0thet(it) c write (iout,*) "thet_pred_mean",thet_pred_mean C Derivatives of the "mean" values in gamma1 and gamma2. - dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss - dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss + dthetg1=(-athet(1,it,ichir1,ichir2)*y(2) + &+athet(2,it,ichir1,ichir2)*y(1))*ss + dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2) + & +bthet(2,it,ichir1,ichir2)*z(1))*ss + if (it.eq.10) then + dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2) + &+athet(2,itype1,ichir11,ichir12)*y(1))*ss + dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2) + & +bthet(2,itype2,ichir21,ichir22)*z(1))*ss + endif if (theta(i).gt.pi-delta) then call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0, & E_tc0) @@ -3250,12 +3672,41 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai +c write (iout,'(a6,i5,0pf7.3,f7.3,i5)') +c & 'ebend',i,ethetai,theta(i),itype(i) c write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i), c & rad2deg*phii,rad2deg*phii1,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 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) - 1215 continue +c 1215 continue + enddo + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +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 enddo C Ufff.... We've done all this!!! return @@ -3370,7 +3821,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation). end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -3390,6 +3841,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@ -3398,37 +3850,53 @@ C etheta=0.0D0 c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) do i=ithet_start,ithet_end - if (itype(i-1).eq.21) cycle +C if (i.eq.2) cycle +C if (itype(i-1).eq.ntyp1) cycle + if (i.le.2) cycle + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle + if (iabs(itype(i+1)).eq.20) iblock=2 + if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp(itype(i-1)) + ityp2=ithetyp((itype(i-1))) do k=1,nntheterm coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.21) then + if (i.eq.3) then + phii=0.0d0 + ityp1=nthetyp+1 + do k=1,nsingle + cosph1(k)=0.0d0 + sinph1(k)=0.0d0 + enddo + else + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 #else phii=phi(i) #endif - ityp1=ithetyp(itype(i-2)) + ityp1=ithetyp((itype(i-2))) do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) enddo else phii=0.0d0 - ityp1=nthetyp+1 +c ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.21) then + endif + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -3436,14 +3904,15 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) #else phii1=phi(i+1) #endif - ityp3=ithetyp(itype(i)) + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) enddo else phii1=0.0d0 - ityp3=nthetyp+1 +c ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -3452,7 +3921,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) c write (iout,*) "i",i," ityp1",itype(i-2),ityp1, c & " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3 c call flush(iout) - ethetai=aa0thet(ityp1,ityp2,ityp3) + ethetai=aa0thet(ityp1,ityp2,ityp3,iblock) do k=1,ndouble do l=1,k-1 ccl=cosph1(l)*cosph2(k-l) @@ -3474,11 +3943,12 @@ c call flush(iout) enddo endif do k=1,ntheterm - ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k) - dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3) + ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k) + dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock) & *coskt(k) if (lprn) - & write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3), + & write (iout,*) "k",k," + & aathet",aathet(k,ityp1,ityp2,ityp3,iblock), & " ethetai",ethetai enddo if (lprn) then @@ -3497,24 +3967,24 @@ c call flush(iout) endif do m=1,ntheterm2 do k=1,nsingle - aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k) - & +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k) - & +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k) - & +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k) + aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k) + & +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k) + & +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k) + & +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k) ethetai=ethetai+sinkt(m)*aux dethetai=dethetai+0.5d0*m*aux*coskt(m) dephii=dephii+k*sinkt(m)*( - & ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)- - & bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)) + & ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)- + & bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)) dephii1=dephii1+k*sinkt(m)*( - & eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)- - & ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k)) + & eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)- + & ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)) if (lprn) & write (iout,*) "m",m," k",k," bbthet", - & bbthet(k,m,ityp1,ityp2,ityp3)," ccthet", - & ccthet(k,m,ityp1,ityp2,ityp3)," ddthet", - & ddthet(k,m,ityp1,ityp2,ityp3)," eethet", - & eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai + & bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet", + & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet", + & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet", + & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai enddo enddo if (lprn) @@ -3522,28 +3992,29 @@ c call flush(iout) do m=1,ntheterm3 do k=2,ndouble do l=1,k-1 - aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+ - & ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+ - & ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l) + aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+ + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+ + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l) ethetai=ethetai+sinkt(m)*aux dethetai=dethetai+0.5d0*m*coskt(m)*aux dephii=dephii+l*sinkt(m)*( - & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)- - & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+ - & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)- + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+ + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)) dephii1=dephii1+(k-l)*sinkt(m)*( - & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+ - & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)- - & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+ + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)- + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)) if (lprn) then write (iout,*) "m",m," k",k," l",l," ffthet", - & ffthet(l,k,m,ityp1,ityp2,ityp3), - & ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet", - & ggthet(l,k,m,ityp1,ityp2,ityp3), - & ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai + & ffthet(l,k,m,ityp1,ityp2,ityp3,iblock), + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet", + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock), + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock), + & " ethetai",ethetai write (iout,*) cosph1ph2(l,k)*sinkt(m), & cosph1ph2(k,l)*sinkt(m), & sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m) @@ -3558,7 +4029,36 @@ c call flush(iout) etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai +c gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai + enddo +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +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) +C endif enddo return end @@ -3586,14 +4086,14 @@ C ALPHA and OMEGA. common /sccalc/ time11,time12,time112,theti,it,nlobit delta=0.02d0*pi escloc=0.0D0 -c write (iout,'(a)') 'ESC' +C write (iout,*) 'ESC' do i=loc_start,loc_end it=itype(i) - if (it.eq.21) cycle + if (it.eq.ntyp1) cycle if (it.eq.10) goto 1 - nlobit=nlob(it) + nlobit=nlob(iabs(it)) c print *,'i=',i,' it=',it,' nlobit=',nlobit -c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad +C write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad theti=theta(i+1)-pipol x(1)=dtan(theti) x(2)=alph(i) @@ -3629,8 +4129,8 @@ c write (iout,*) "i",i," x",x(1),x(2),x(3) dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) enddo dersc(2)=dersc(2)+ssd*(escloci-esclocbi) -c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, -c & esclocbi,ss,ssd + write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, + & esclocbi,ss,ssd escloci=ss*escloci+(1.0d0-ss)*esclocbi c escloci=esclocbi c write (iout,*) escloci @@ -3664,15 +4164,17 @@ c write (iout,*) escloci enddo dersc(2)=dersc(2)+ssd*(escloci-esclocbi) c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, -c & esclocbi,ss,ssd +c & esclocbi,ss,ssd escloci=ss*escloci+(1.0d0-ss)*esclocbi -c write (iout,*) escloci +C write (iout,*) 'i=',i, escloci else call enesc(x,escloci,dersc,ddummy,.false.) endif escloc=escloc+escloci -c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc +C write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc + write (iout,'(a6,i5,0pf7.3)') + & 'escloc',i,escloci gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & wscloc*dersc(1) @@ -3746,7 +4248,7 @@ C Compute the contribution to SC energy and derivatives do iii=-1,1 do j=1,nlobit - expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin) cd print *,'j=',j,' expfac=',expfac escloc_i=escloc_i+expfac do k=1,3 @@ -3827,7 +4329,7 @@ C Compute the contribution to SC energy and derivatives dersc12=0.0d0 do j=1,nlobit - expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin) escloc_i=escloc_i+expfac do k=1,2 dersc(k)=dersc(k)+Ax(k,j)*expfac @@ -3882,7 +4384,7 @@ C delta=0.02d0*pi escloc=0.0D0 do i=loc_start,loc_end - if (itype(i).eq.21) cycle + if (itype(i).eq.ntyp1) cycle costtab(i+1) =dcos(theta(i+1)) sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1)) cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1))) @@ -3891,7 +4393,7 @@ C cosfac=dsqrt(cosfac2) sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) - it=itype(i) + it=iabs(itype(i)) if (it.eq.10) goto 1 c C Compute the axes of tghe local cartesian coordinates system; store in @@ -3909,7 +4411,7 @@ C & dc_norm(3,i+nres) y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac enddo do j = 1,3 - z_prime(j) = -uz(j,i-1) + z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i))) enddo c write (2,*) "i",i c write (2,*) "x_prime",(x_prime(j),j=1,3) @@ -3941,7 +4443,7 @@ C C Compute the energy of the ith side cbain C c write (2,*) "xx",xx," yy",yy," zz",zz - it=itype(i) + it=iabs(itype(i)) do j = 1,65 x(j) = sc_parmin(j,it) enddo @@ -3949,7 +4451,7 @@ c write (2,*) "xx",xx," yy",yy," zz",zz Cc diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0d0,itype(i))*dsin(alph(2))*dsin(omeg(2)) write(2,'(3f8.1,3f9.3,1x,3f9.3)') & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz, & xx1,yy1,zz1 @@ -3992,6 +4494,8 @@ c & dscp1,dscp2,sumene c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene c write (2,*) "escloc",escloc +c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i), +c & zz,xx,yy if (.not. calc_grad) goto 1 #ifdef DEBUG C @@ -4120,8 +4624,10 @@ c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i) dZZ_Ci1(k)=0.0d0 dZZ_Ci(k)=0.0d0 do j=1,3 - dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres) - dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres) + dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) + & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) + dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) + & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) enddo dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres)) @@ -4273,8 +4779,8 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.21 .or. itype(i-1).eq.21 - & .or. itype(i).eq.21) cycle + if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1) cycle itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) phii=phi(i) @@ -4322,15 +4828,16 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) difi=phii-phi0(i) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 endif -! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, -! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) +C write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih", +C & i,itori,rad2deg*phii, +C & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) enddo ! write (iout,*) 'edihcnstr',edihcnstr return @@ -4358,17 +4865,25 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.21 .or. itype(i-1).eq.21 - & .or. itype(i).eq.21) cycle + if (i.le.2) cycle + if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle +C if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 +C & .or. itype(i).eq.ntyp1) cycle if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215 + if (iabs(itype(i)).eq.20) then + iblock=2 + else + iblock=1 + endif itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 C Regular cosine and sine terms - do j=1,nterm(itori,itori1) - v1ij=v1(j,itori,itori1) - v2ij=v2(j,itori,itori1) + do j=1,nterm(itori,itori1,iblock) + v1ij=v1(j,itori,itori1,iblock) + v2ij=v2(j,itori,itori1,iblock) cosphi=dcos(j*phii) sinphi=dsin(j*phii) etors=etors+v1ij*cosphi+v2ij*sinphi @@ -4381,22 +4896,24 @@ C [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1 C cosphi=dcos(0.5d0*phii) sinphi=dsin(0.5d0*phii) - do j=1,nlor(itori,itori1) + do j=1,nlor(itori,itori1,iblock) vl1ij=vlor1(j,itori,itori1) vl2ij=vlor2(j,itori,itori1) vl3ij=vlor3(j,itori,itori1) pom=vl2ij*cosphi+vl3ij*sinphi pom1=1.0d0/(pom*pom+1.0d0) etors=etors+vl1ij*pom1 +c if (energy_dec) etors_ii=etors_ii+ +c & vl1ij*pom1 pom=-pom*pom1*pom1 gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom enddo C Subtract the constant term - etors=etors-v0(itori,itori1) + etors=etors-v0(itori,itori1,iblock) if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, - & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) + & (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) 1215 continue @@ -4410,21 +4927,24 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) edihi=0.0d0 if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 - edihi=0.25d0*ftors*difi**4 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 + edihi=0.25d0*ftors(i)*difi**4 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 - edihi=0.25d0*ftors*difi**4 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 + edihi=0.25d0*ftors(i)*difi**4 else difi=0.0d0 endif + write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih", + & i,itori,rad2deg*phii, + & rad2deg*difi,0.25d0*ftors(i)*difi**4 c write (iout,'(2i5,4f10.5,e15.5)') i,itori,phii,phi0(i),difi, c & drange(i),edihi ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, -! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) +! & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) enddo ! write (iout,*) 'edihcnstr',edihcnstr return @@ -4452,8 +4972,12 @@ C Set lprn=.true. for debugging c lprn=.true. etors_d=0.0D0 do i=iphi_start,iphi_end-1 - if (itype(i-2).eq.21 .or. itype(i-1).eq.21 - & .or. itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (i.le.3) cycle +C if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 +C & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or. + & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or. + & (itype(i+1).eq.ntyp1)) cycle if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0) & goto 1215 itori=itortyp(itype(i-2)) @@ -4463,12 +4987,14 @@ c lprn=.true. phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 + iblock=1 + if (iabs(itype(i+1)).eq.20) iblock=2 C Regular cosine and sine terms - do j=1,ntermd_1(itori,itori1,itori2) - v1cij=v1c(1,j,itori,itori1,itori2) - v1sij=v1s(1,j,itori,itori1,itori2) - v2cij=v1c(2,j,itori,itori1,itori2) - v2sij=v1s(2,j,itori,itori1,itori2) + do j=1,ntermd_1(itori,itori1,itori2,iblock) + v1cij=v1c(1,j,itori,itori1,itori2,iblock) + v1sij=v1s(1,j,itori,itori1,itori2,iblock) + v2cij=v1c(2,j,itori,itori1,itori2,iblock) + v2sij=v1s(2,j,itori,itori1,itori2,iblock) cosphi1=dcos(j*phii) sinphi1=dsin(j*phii) cosphi2=dcos(j*phii1) @@ -4478,12 +5004,12 @@ C Regular cosine and sine terms gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo - do k=2,ntermd_2(itori,itori1,itori2) + do k=2,ntermd_2(itori,itori1,itori2,iblock) do l=1,k-1 - v1cdij = v2c(k,l,itori,itori1,itori2) - v2cdij = v2c(l,k,itori,itori1,itori2) - v1sdij = v2s(k,l,itori,itori1,itori2) - v2sdij = v2s(l,k,itori,itori1,itori2) + v1cdij = v2c(k,l,itori,itori1,itori2,iblock) + v2cdij = v2c(l,k,itori,itori1,itori2,iblock) + v1sdij = v2s(k,l,itori,itori1,itori2,iblock) + v2sdij = v2s(l,k,itori,itori1,itori2,iblock) cosphi1p2=dcos(l*phii+(k-l)*phii1) cosphi1m2=dcos(l*phii-(k-l)*phii1) sinphi1p2=dsin(l*phii+(k-l)*phii1) @@ -4493,7 +5019,7 @@ C Regular cosine and sine terms gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 & -v1cdij*sinphi1p2-v2cdij*sinphi1m2) gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 - & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) + & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) enddo enddo gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*fact2*gloci1 @@ -4532,26 +5058,50 @@ C Set lprn=.true. for debugging c lprn=.true. c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor esccor=0.0D0 - do i=iphi_start,iphi_end - if (itype(i-2).eq.21 .or. itype(i-1).eq.21) cycle + do i=itau_start,itau_end + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle esccor_ii=0.0D0 - itori=itype(i-2) - itori1=itype(i-1) + isccori=isccortyp(itype(i-2)) + isccori1=isccortyp(itype(i-1)) phii=phi(i) + do intertyp=1,3 !intertyp +cc Added 09 May 2012 (Adasko) +cc Intertyp means interaction type of backbone mainchain correlation: +c 1 = SC...Ca...Ca...Ca +c 2 = Ca...Ca...Ca...SC +c 3 = SC...Ca...Ca...SCi gloci=0.0D0 - do j=1,nterm_sccor - v1ij=v1sccor(j,itori,itori1) - v2ij=v2sccor(j,itori,itori1) - cosphi=dcos(j*phii) - sinphi=dsin(j*phii) - esccor=esccor+v1ij*cosphi+v2ij*sinphi - gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) - enddo + if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. + & (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-1).eq.ntyp1))) + & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) + & .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1) + & .or.(itype(i).eq.ntyp1))) + & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. + & (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-3).eq.ntyp1)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) + & cycle + do j=1,nterm_sccor(isccori,isccori1) + v1ij=v1sccor(j,intertyp,isccori,isccori1) + v2ij=v2sccor(j,intertyp,isccori,isccori1) + cosphi=dcos(j*tauangle(intertyp,i)) + sinphi=dsin(j*tauangle(intertyp,i)) + esccor=esccor+v1ij*cosphi+v2ij*sinphi + gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) + enddo +C write (iout,*)"EBACK_SC_COR",esccor,i +c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp, +c & nterm_sccor(isccori,isccori1),isccori,isccori1 +c gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, - & (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6) - gsccor_loc(i-3)=gloci + & (v1sccor(j,1,itori,itori1),j=1,6) + & ,(v2sccor(j,1,itori,itori1),j=1,6) +c gsccor_loc(i-3)=gloci + enddo !intertyp enddo return end @@ -4660,9 +5210,9 @@ c------------------------------------------------------------------------------ integer dimen1,dimen2,atom,indx double precision buffer(dimen1,dimen2) double precision zapas - common /contacts_hb/ zapas(3,20,maxres,7), - & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres), - & num_cont_hb(maxres),jcont_hb(20,maxres) + common /contacts_hb/ zapas(3,ntyp,maxres,7), + & facont_hb(ntyp,maxres),ees0p(ntyp,maxres),ees0m(ntyp,maxres), + & num_cont_hb(maxres),jcont_hb(ntyp,maxres) num_kont=num_cont_hb(atom) do i=1,num_kont do k=1,7 @@ -4685,9 +5235,10 @@ c------------------------------------------------------------------------------ integer dimen1,dimen2,atom,indx double precision buffer(dimen1,dimen2) double precision zapas - common /contacts_hb/ zapas(3,20,maxres,7), - & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres), - & num_cont_hb(maxres),jcont_hb(20,maxres) + common /contacts_hb/ zapas(3,ntyp,maxres,7), + & facont_hb(ntyp,maxres),ees0p(ntyp,maxres), + & ees0m(ntyp,maxres), + & num_cont_hb(maxres),jcont_hb(ntyp,maxres) num_kont=buffer(1,indx+26) num_kont_old=num_cont_hb(atom) num_cont_hb(atom)=num_kont+num_kont_old @@ -6426,7 +6977,7 @@ c---------------------------------------------------------------------------- include 'COMMON.GEO' logical swap double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), - & auxvec1(2),auxvec2(1),auxmat1(2,2) + & auxvec1(2),auxvec2(2),auxmat1(2,2) logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC @@ -7303,6 +7854,125 @@ cd write (2,*) 'eel_turn6',ekont*eel_turn6 return end crc------------------------------------------------- +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + subroutine Eliptransfer(eliptran) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' +C this is done by Adasko +C print *,"wchodze" +C structure of box: +C water +C--bordliptop-- buffore starts +C--bufliptop--- here true lipid starts +C lipid +C--buflipbot--- lipid ends buffore starts +C--bordlipbot--buffore ends + eliptran=0.0 + do i=1,nres +C do i=1,1 + if (itype(i).eq.ntyp1) cycle + + positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize)) + if (positi.le.0) positi=positi+boxzsize +C print *,i +C first for peptide groups +c for each residue check if it is in lipid or lipid water border area + if ((positi.gt.bordlipbot) + &.and.(positi.lt.bordliptop)) then +C the energy transfer exist + if (positi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 +C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 +C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran +C print *, "doing sscalefor top part" +C print *,i,sslip,fracinbuf,ssgradlip + else + eliptran=eliptran+pepliptran +C print *,"I am in true lipid" + endif +C else +C eliptran=elpitran+0.0 ! I am in water + endif + enddo +C print *, "nic nie bylo w lipidzie?" +C now multiply all by the peptide group transfer factor +C eliptran=eliptran*pepliptran +C now the same for side chains +CV do i=1,1 + do i=1,nres + if (itype(i).eq.ntyp1) cycle + positi=(mod(c(3,i+nres),boxzsize)) + if (positi.le.0) positi=positi+boxzsize +C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop +c for each residue check if it is in lipid or lipid water border area +C respos=mod(c(3,i+nres),boxzsize) +C print *,positi,bordlipbot,buflipbot + if ((positi.gt.bordlipbot) + & .and.(positi.lt.bordliptop)) then +C the energy transfer exist + if (positi.lt.buflipbot) then + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0- + &((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + eliptran=eliptran+liptranene(itype(i)) +C print *,"I am in true lipid" + endif + endif ! if in lipid or buffor +C else +C eliptran=elpitran+0.0 ! I am in water + enddo + return + end + + +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + SUBROUTINE MATVEC2(A1,V1,V2) implicit real*8 (a-h,o-z) include 'DIMENSIONS' @@ -7436,4 +8106,64 @@ C----------------------------------------------------------------------------- scalar=sc return end +C----------------------------------------------------------------------- + double precision function sscale(r) + double precision r,gamm + include "COMMON.SPLITELE" + if(r.lt.r_cut-rlamb) then + sscale=1.0d0 + else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then + gamm=(r-(r_cut-rlamb))/rlamb + sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0) + else + sscale=0d0 + endif + return + end +C----------------------------------------------------------------------- +C----------------------------------------------------------------------- + double precision function sscagrad(r) + double precision r,gamm + include "COMMON.SPLITELE" + if(r.lt.r_cut-rlamb) then + sscagrad=0.0d0 + else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then + gamm=(r-(r_cut-rlamb))/rlamb + sscagrad=gamm*(6*gamm-6.0d0)/rlamb + else + sscagrad=0.0d0 + endif + return + end +C----------------------------------------------------------------------- +C----------------------------------------------------------------------- + double precision function sscalelip(r) + double precision r,gamm + include "COMMON.SPLITELE" +C if(r.lt.r_cut-rlamb) then +C sscale=1.0d0 +C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then +C gamm=(r-(r_cut-rlamb))/rlamb + sscalelip=1.0d0+r*r*(2*r-3.0d0) +C else +C sscale=0d0 +C endif + return + end +C----------------------------------------------------------------------- + double precision function sscagradlip(r) + double precision r,gamm + include "COMMON.SPLITELE" +C if(r.lt.r_cut-rlamb) then +C sscagrad=0.0d0 +C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then +C gamm=(r-(r_cut-rlamb))/rlamb + sscagradlip=r*(6*r-6.0d0) +C else +C sscagrad=0.0d0 +C endif + return + end + +C-----------------------------------------------------------------------