X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fenergy_p_new-sep_barrier.F;h=f92aebb9ad0031a150c3f194a461640a8eb032fb;hb=4627052ad84994541a6f518b197289553d91edc0;hp=96f77770fdf54e59677d209122ae088cca018393;hpb=566219e74d5310efba22c5a154c9075d971f6603;p=unres.git diff --git a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F index 96f7777..f92aebb 100644 --- a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F @@ -84,14 +84,19 @@ c include 'COMMON.CONTACTS' integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sss1,sssgrad1,sqrij - double precision sscale,sscagrad + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip + double precision sscale,sscagrad,sscagradlip,sscalelip double precision boxshift + double precision gg_lipi(3),gg_lipj(3) c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -99,6 +104,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C C Calculate SC interaction energy. C @@ -112,6 +118,11 @@ c do j=istart(i,iint),iend(i,iint) 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -127,6 +138,7 @@ c do j=istart(i,iint),iend(i,iint) if (sss.lt.1.0d0) then rrij=1.0D0/rij fac=rrij**expon2 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=e1+e2 @@ -140,11 +152,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo endif c enddo ! j @@ -192,14 +209,19 @@ c include 'COMMON.CONTACTS' integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift + double precision gg_lipi(3),gg_lipj(3) c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -207,6 +229,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C Change 12/1/95 num_conti=0 C @@ -222,6 +245,11 @@ c do j=istart(i,iint),iend(i,iint) 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -235,6 +263,7 @@ C Change 12/1/95 to calculate four-body interactions rrij=1.0D0/rij eps0ij=eps(itypi,itypj) fac=rrij**expon2 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=e1+e2 @@ -246,11 +275,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo endif c enddo ! j @@ -296,14 +330,19 @@ C double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift + double precision gg_lipi(3),gg_lipj(3) c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -311,6 +350,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C C Calculate SC interaction energy. C @@ -322,6 +362,11 @@ c do j=istart(i,iint),iend(i,iint) 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -339,6 +384,7 @@ c do j=istart(i,iint),iend(i,iint) & sscagrad(rij/sigma(itypi,itypj),r_cut_respa) r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=e_augm+e1+e2 @@ -360,11 +406,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo endif c enddo ! j @@ -401,14 +452,19 @@ C double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift + double precision gg_lipi(3),gg_lipj(3) c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -416,6 +472,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C C Calculate SC interaction energy. C @@ -427,6 +484,11 @@ c do j=istart(i,iint),iend(i,iint) 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -439,6 +501,7 @@ c do j=istart(i,iint),iend(i,iint) if (sss.gt.0.0d0) then r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=e_augm+e1+e2 @@ -459,11 +522,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo endif c enddo ! j @@ -501,13 +569,16 @@ C integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision sss1,sssgrad1 - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift c double precision rrsave(maxdim) logical lprn evdw=0.0D0 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c if (icall.eq.0) then c lprn=.true. c else @@ -515,9 +586,9 @@ c else c endif ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -525,6 +596,7 @@ c do i=iatsc_s,iatsc_e 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) @@ -553,6 +625,11 @@ c dscj_inv=dsc_inv(itypj) 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -573,6 +650,7 @@ 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 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -600,6 +678,12 @@ C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad_scale((1.0d0-sss)*sss1) @@ -633,13 +717,16 @@ C double precision evdw integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift c double precision rrsave(maxdim) logical lprn evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 c if (icall.eq.0) then c lprn=.true. c else @@ -647,9 +734,9 @@ c else c endif ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -657,6 +744,7 @@ c do i=iatsc_s,iatsc_e 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) @@ -685,6 +773,11 @@ c dscj_inv=dsc_inv(itypj) 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -702,6 +795,7 @@ 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 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -728,6 +822,17 @@ C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi + do k=1,3 + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) + enddo C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad_scale(sss) @@ -767,18 +872,22 @@ C & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision subchap,sss1,sssgrad1 - double precision boxshift + double precision boxshift,rij1 evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + if (energy_dec) + & write(2,*) "g_listscsc_start_long,g_listscsc_end_long", + & g_listscsc_start_long,g_listscsc_end_long + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -823,9 +932,9 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(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 + & +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 + & +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) @@ -834,9 +943,13 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss1=sscale(1.0d0/rij,r_cut_int) + rij1=1.0d0/rij +c sss1=sscale(1.0d0/rij,r_cut_int) + sss1=sscale(rij1,r_cut_int) if (sss1.eq.0.0d0) cycle - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) + rij1=rij1/sigmaii(itypi,itypj) + sss=sscale(rij1,r_cut_respa) +c sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) if (sss.lt.1.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. @@ -861,6 +974,7 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -881,8 +995,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a6,2i5,4f10.5)') - & 'evdw',i,j,rij,sss,sss1,evdwij + if (energy_dec) write (iout,'(a,2i5,5f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss1,sss,sslipi,sslipj,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -895,8 +1009,13 @@ 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 + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi + C Calculate angular part of the gradient. call sc_grad_scale((1.0d0-sss)*sss1) endif @@ -914,6 +1033,7 @@ C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne potential of interaction. C implicit none + include 'mpif.h' include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -926,6 +1046,7 @@ C include 'COMMON.CALC' include 'COMMON.CONTROL' include "COMMON.SPLITELE" + include 'COMMON.TIME1' logical lprn double precision evdw integer itypi,itypj,itypi1,iint,ind,ikont @@ -934,17 +1055,23 @@ C & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision boxshift + double precision time01 +c time01=MPI_Wtime() evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0D0 + gg_lipj=0.0d0 lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + if (energy_dec) + & write(2,*) "g_listscsc_start_short,g_listscsc_end_short", + & g_listscsc_start_short,g_listscsc_end_short + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -989,9 +1116,14 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(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 + & +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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 +c write (iout,*) "aa bb",aa_lip(itypi,itypj), +c & bb_lip(itypi,itypj),aa_aq(itypi,itypj), +c & bb_aq(itypi,itypj),aa,bb +c write (iout,*) (sslipi+sslipj)/2.0d0, +c & (2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1001,8 +1133,8 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) - sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) if (sss.gt.0.0d0) then + sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) C Calculate angle-dependent terms of energy and contributions to their C derivatives. @@ -1024,6 +1156,7 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1044,8 +1177,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw',i,j,evdwij + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1057,14 +1190,20 @@ 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 + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi +c write (iout,*) "gglip",i,j,gg_lipi,gg_lipj C Calculate angular part of the gradient. call sc_grad_scale(sss) endif c enddo ! j c enddo ! iint enddo ! i +c time_evdw_short=time_evdw_short+MPI_Wtime()-time01 c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return @@ -1099,15 +1238,16 @@ C double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision sss1,sssgrad1 evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 lprn=.false. c if (icall.eq.0) lprn=.true. ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1183,6 +1323,7 @@ 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 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1215,6 +1356,12 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate angular part of the gradient. call sc_grad_scale((1.0d0-sss)*sss1) endif @@ -1253,15 +1400,16 @@ C double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision boxshift evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 lprn=.false. c if (icall.eq.0) lprn=.true. ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1333,6 +1481,7 @@ 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 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1363,6 +1512,12 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate angular part of the gradient. call sc_grad_scale(sss) endif @@ -1417,6 +1572,7 @@ c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv C C Calculate the components of the gradient in DC and X C +c write (iout,*) "scgrad gglip",i,j,gg_lipi,gg_lipj do l=1,3 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) @@ -1465,6 +1621,8 @@ C common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -1559,6 +1717,7 @@ C & .or. itype(i+4).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 call eelecij_scale(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -1583,6 +1742,7 @@ C & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -1597,6 +1757,9 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c c do i=iatel_s,iatel_e + if (energy_dec) + & write(iout,*) "g_listpp_start,g_listpp_end", + & g_listpp_start,g_listpp_end do ikont=g_listpp_start,g_listpp_end i=newcontlistppi(ikont) j=newcontlistppj(ikont) @@ -1614,6 +1777,7 @@ C & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend #ifdef FOURBODY num_conti=num_cont_hb(i) @@ -1696,7 +1860,9 @@ C------------------------------------------------------------------------------- double precision sscale,sscagrad double precision scalar double precision boxshift - + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij, + & faclipij2 + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -1729,6 +1895,9 @@ C print *,"WCHODZE2" yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0 + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -1766,17 +1935,18 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions eesij=el1+el2 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*sss1 - evdw1=evdw1+evdwij*(1.0d0-sss)*sss1 + ees=ees+eesij*sss1*faclipij2 + evdw1=evdw1+evdwij*(1.0d0-sss)*sss1*faclipij2 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 if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3,2f7.3)') - & 'evdw1',i,j,evdwij,sss,sss1 - write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij + write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,5f10.5)') + & 'evdw1',i,j,evdwij,iteli,itelj,aaa,sss,sss1,sssgrad,sssgrad1,rij + write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij,faclipij2 endif C @@ -1794,7 +1964,8 @@ c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) * * Radial derivatives. First process both termini of the fragment (i,j) * - aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij +c old aux=(facel+sssgrad1*(1.0d0-sss)*eesij*rmij)*faclipij2 + aux=(facel+sssgrad1*eesij*rmij)*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) ggg(1)=aux*xj ggg(2)=aux*yj @@ -1867,6 +2038,10 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo + gelc_long(3,j)=gelc_long(3,j)+ + & ssgradlipj*eesij/2.0d0*lipscale**2*sss1 + gelc_long(3,i)=gelc_long(3,i)+ + & ssgradlipi*eesij/2.0d0*lipscale**2*sss1 c gelc_long(3,i)=gelc_long(3,i)+ c ssgradlipi*eesij/2.0d0*lipscale**2*sss1 * @@ -1877,9 +2052,9 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - facvdw=facvdw+ - & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij -c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + facvdw=(facvdw+ + &(-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij) + & *faclipij2 ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj @@ -1893,6 +2068,11 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +!C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss1*(1.0d0-sss)*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss1*(1.0d0-sss)*ssgradlipi*evdwij/2.0d0*lipscale**2 * * Loop over residues i+1 thru j-1. * @@ -1917,8 +2097,8 @@ c facel=el1+eesij * * Radial derivatives. First process both termini of the fragment (i,j) * - aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj)) - & *eesij*rmij + aux=(fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj)) + & *eesij*rmij)*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) ggg(1)=aux*xj ggg(2)=aux*yj @@ -1957,6 +2137,10 @@ C ggg(3)=facvdw*zj gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss1*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss1*ssgradlipi*evdwij/2.0d0*lipscale**2 #endif * * Angular part @@ -1974,7 +2158,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1 - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) enddo @@ -1996,13 +2180,13 @@ cgrad enddo gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1 - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1 - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) @@ -2216,7 +2400,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)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij eel_loc=eel_loc+eel_loc_ij if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. @@ -2268,7 +2452,7 @@ C & *2.0 & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -2284,7 +2468,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -2297,7 +2481,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij geel_loc_ji= & +a22*gmuji2(1) @@ -2309,14 +2493,14 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij #endif cC Paral derivatives in virtual-bond dihedral angles gamma if (i.gt.1) & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) @@ -2324,7 +2508,7 @@ c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) @@ -2336,7 +2520,7 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) @@ -2346,6 +2530,10 @@ cgrad ghalf=0.5d0*ggg(l) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo + gel_loc_long(3,j)=gel_loc_long(3,j)+ + & ssgradlipj*eel_loc_ij/2.0d0*lipscale/faclipij + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -2363,22 +2551,22 @@ c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut do l=1,3 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo @@ -2632,17 +2820,22 @@ c write (iout,*) "evdwpp_short" double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, & dist_temp, dist_init,sss_grad double precision sscale,sscagrad + double precision sslipi,ssgradlipi,sslipj,ssgradlipj double precision boxshift integer ikont + double precision faclipij2 evdw1=0.0D0 C print *,"WCHODZE" c write (iout,*) "iatel_s_vdw",iatel_s_vdw, c & " iatel_e_vdw",iatel_e_vdw c call flush(iout) c do i=iatel_s_vdw,iatel_e_vdw - do ikont=g_listpp_vdw_start,g_listpp_vdw_end - i=newcontlistpp_vdwi(ikont) - j=newcontlistpp_vdwj(ikont) + if (energy_dec) + & write(iout,*) "g_listpp_vdw_start_short,g_listpp_vdw_end_short", + & g_listpp_vdw_start_short,g_listpp_vdw_end_short + do ikont=g_listpp_vdw_start_short,g_listpp_vdw_end_short + i=newcontlistpp_vdwi_short(ikont) + j=newcontlistpp_vdwj_short(ikont) if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -2654,6 +2847,7 @@ c do i=iatel_s_vdw,iatel_e_vdw ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i), c & ' ielend',ielend_vdw(i) @@ -2676,6 +2870,8 @@ c do j=ielstart_vdw(i),ielend_vdw(i) yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -2698,16 +2894,17 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions if (energy_dec) then write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss endif - evdw1=evdw1+evdwij*sss + evdw1=evdw1+evdwij*sss*faclipij2 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)') & 'evdw1_sum',i,j,evdw1 C C Calculate contributions to the Cartesian gradient. C - facvdw=-6*rrmij*(ev1+evdwij)*sss - ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) - ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) - ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) + facvdw=(-6*rrmij*(ev1+evdwij)*sss+sssgrad*rmij*evdwij/ + & rpp(iteli,itelj))*faclipij2 + 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 @@ -2715,6 +2912,11 @@ C ggg(3)=facvdw*zj gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +!C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 endif c enddo ! j enddo ! i @@ -2758,9 +2960,12 @@ c if (lprint_short) c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s, c & ' iatscp_e=',iatscp_e c do i=iatscp_s,iatscp_e - do ikont=g_listscp_start,g_listscp_end - i=newcontlistscpi(ikont) - j=newcontlistscpj(ikont) + if (energy_dec) + & write(iout,*)"g_listscp_start_long,g_listscp_end_long", + & g_listscp_start_long,g_listscp_end_long + do ikont=g_listscp_start_long,g_listscp_end_long + i=newcontlistscpi_long(ikont) + j=newcontlistscpj_long(ikont) if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) @@ -2884,14 +3089,20 @@ C double precision ggg(3) double precision sscale,sscagrad double precision boxshift + integer ikont evdw2=0.0D0 evdw2_14=0.0d0 cd print '(a)','Enter ESCP' c if (lprint_short) c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s, c & ' iatscp_e=',iatscp_e - if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb - do i=iatscp_s,iatscp_e +c if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb + if (energy_dec) + & write(iout,*) "g_listscp_start_short,g_listscp_end_short", + & g_listscp_start_short,g_listscp_end_short + do ikont=g_listscp_start_short,g_listscp_end_short + i=newcontlistscpi_short(ikont) + j=newcontlistscpj_short(ikont) if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) @@ -2902,9 +3113,9 @@ c & ' iatscp_e=',iatscp_e c if (lprint_short) c & write (iout,*) "i",i," itype",itype(i),itype(i+1), c & " nscp_gr",nscp_gr(i) - do iint=1,nscp_gr(i) - - do j=iscpstart(i,iint),iscpend(i,iint) +c do iint=1,nscp_gr(i) +c +c do j=iscpstart(i,iint),iscpend(i,iint) itypj=iabs(itype(j)) c if (lprint_short) c & write (iout,*) "j",j," itypj",itypj @@ -2968,9 +3179,9 @@ c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo endif - enddo +c enddo - enddo ! iint +c enddo ! iint enddo ! i do i=1,nct do j=1,3