X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-SAXS%2Fenergy_p_new-sep_barrier.F;h=261c06c8ebb92da811ff965174cc4848af29ece8;hb=1346fb3959c2eb0a370b11bc6ccad5e4cca27ec9;hp=1f00b2be0cab7209eb77c47484020dfa2f69ca5a;hpb=5836ecdab5a8b95f079bbf6e07374dee3fce8a26;p=unres.git diff --git a/source/unres/src_MD-M-SAXS/energy_p_new-sep_barrier.F b/source/unres/src_MD-M-SAXS/energy_p_new-sep_barrier.F index 1f00b2b..261c06c 100644 --- a/source/unres/src_MD-M-SAXS/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M-SAXS/energy_p_new-sep_barrier.F @@ -1,5 +1,6 @@ C----------------------------------------------------------------------- double precision function sscalelip(r) + implicit none double precision r,gamm include "COMMON.SPLITELE" C if(r.lt.r_cut-rlamb) then @@ -14,6 +15,7 @@ C endif end C----------------------------------------------------------------------- double precision function sscagradlip(r) + implicit none double precision r,gamm include "COMMON.SPLITELE" C if(r.lt.r_cut-rlamb) then @@ -28,8 +30,9 @@ C endif end C----------------------------------------------------------------------- - double precision function sscale(r) - double precision r,gamm + double precision function sscale(r,r_cut) + implicit none + double precision r,r_cut,gamm include "COMMON.SPLITELE" if(r.lt.r_cut-rlamb) then sscale=1.0d0 @@ -42,9 +45,9 @@ C----------------------------------------------------------------------- return end C----------------------------------------------------------------------- -C----------------------------------------------------------------------- - double precision function sscagrad(r) - double precision r,gamm + double precision function sscagrad(r,r_cut) + implicit none + double precision r,r_cut,gamm include "COMMON.SPLITELE" if(r.lt.r_cut-rlamb) then sscagrad=0.0d0 @@ -62,9 +65,8 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJ potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' - parameter (accur=1.0d-10) include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' @@ -75,14 +77,20 @@ C include 'COMMON.SBRIDGE' include 'COMMON.NAMES' include 'COMMON.IOUNITS' - include 'COMMON.CONTACTS' - dimension gg(3) + include "COMMON.SPLITELE" +c include 'COMMON.CONTACTS' + double precision gg(3) + double precision evdw,evdwij + integer i,j,k,itypi,itypj,itypi1,num_conti,iint + double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, + & sigij,r0ij,rcut,sss1,sssgrad1,sqrij + double precision sscale,sscagrad c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -93,25 +101,33 @@ 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) + 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 rij=xj*xj+yj*yj+zj*zj - sss=sscale(dsqrt(rij)/sigma(itypi,itypj)) + sqrij=dsqrt(rrij) + eps0ij=eps(itypi,itypj) + sss1=sscale(sqrij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(sqrij,r_cut_int) + sssgrad= + & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa) + sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa) if (sss.lt.1.0d0) then rrij=1.0D0/rij - eps0ij=eps(itypi,itypj) fac=rrij**expon2 e1=fac*fac*aa e2=fac*bb evdwij=e1+e2 - evdw=evdw+(1.0d0-sss)*evdwij + evdw=evdw+(1.0d0-sss)*sss1*evdwij C C Calculate the components of the gradient in DC and X C - fac=-rrij*(e1+evdwij)*(1.0d0-sss) + fac=-rrij*(e1+evdwij)*(1.0d0-sss)*sss1 + & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj) + & +(1.0d0-sss)*sssgrad1)/sqrij gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -148,9 +164,8 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJ potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' - parameter (accur=1.0d-10) include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' @@ -161,14 +176,20 @@ C include 'COMMON.SBRIDGE' include 'COMMON.NAMES' include 'COMMON.IOUNITS' - include 'COMMON.CONTACTS' - dimension gg(3) + include "COMMON.SPLITELE" +c include 'COMMON.CONTACTS' + double precision gg(3) + double precision evdw,evdwij + integer i,j,k,itypi,itypj,itypi1,num_conti,iint + double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, + & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 + double precision sscale,sscagrad c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -181,15 +202,18 @@ 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) + 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 C Change 12/1/95 to calculate four-body interactions rij=xj*xj+yj*yj+zj*zj - sss=sscale(dsqrt(rij)/sigma(itypi,itypj)) + sqrij=dsqrt(rij) + sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa) if (sss.gt.0.0d0) then + sssgrad= + & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa) rrij=1.0D0/rij eps0ij=eps(itypi,itypj) fac=rrij**expon2 @@ -200,7 +224,7 @@ C Change 12/1/95 to calculate four-body interactions C C Calculate the components of the gradient in DC and X C - fac=-rrij*(e1+evdwij)*sss + fac=-rrij*(e1+evdwij)*sss+evdwij*sssgrad/sqrij gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -237,7 +261,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJK potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -247,14 +271,20 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' - dimension gg(3) + include "COMMON.SPLITELE" + double precision gg(3) + double precision evdw,evdwij + integer i,j,k,itypi,itypj,itypi1,iint + 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 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -263,7 +293,7 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -273,8 +303,13 @@ C e_augm=augm(itypi,itypj)*fac_augm r_inv_ij=dsqrt(rrij) rij=1.0D0/r_inv_ij - sss=sscale(rij/sigma(itypi,itypj)) + sss1=sscale(rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(rij,r_cut_int) + sss=sscale(rij/sigma(itypi,itypj),r_cut_respa) if (sss.lt.1.0d0) then + sssgrad= + & 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 e1=fac*fac*aa @@ -287,12 +322,14 @@ 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) - evdw=evdw+(1.0d0-sss)*evdwij + evdw=evdw+(1.0d0-sss)*sss1*evdwij C C Calculate the components of the gradient in DC and X C fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) - fac=fac*(1.0d0-sss) + fac=fac*(1.0d0-sss)*sss1 + & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj) + & +(1.0d0-sss)*sssgrad1)*r_inv_ij gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -330,14 +367,20 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' - dimension gg(3) + include "COMMON.SPLITELE" + double precision gg(3) + double precision evdw,evdwij + integer i,j,k,itypi,itypj,itypi1,iint + 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 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -346,7 +389,7 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -356,7 +399,7 @@ C e_augm=augm(itypi,itypj)*fac_augm r_inv_ij=dsqrt(rrij) rij=1.0D0/r_inv_ij - sss=sscale(rij/sigma(itypi,itypj)) + sss=sscale(rij/sigma(itypi,itypj),r_cut_respa) if (sss.gt.0.0d0) then r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon @@ -375,6 +418,7 @@ C C Calculate the components of the gradient in DC and X C fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) + & +evdwij*sssgrad/sigma(itypi,itypj)*r_inv_ij fac=fac*sss gg(1)=xj*fac gg(2)=yj*fac @@ -403,7 +447,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Berne-Pechukas potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -414,7 +458,14 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include "COMMON.SPLITELE" + integer icall common /srutu/ icall + double precision evdw + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi + double precision sss1,sssgrad1 + double precision sscale,sscagrad c double precision rrsave(maxdim) logical lprn evdw=0.0D0 @@ -427,9 +478,9 @@ c else c endif ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -444,7 +495,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -465,10 +516,13 @@ c dscj_inv=dsc_inv(itypj) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) - + sss1=sscale(1.0d0/rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) if (sss.lt.1.0d0) then - + sssgrad= + & sscagrad(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) + sssgrad1=sscagrad(1.0d0/rij,r_cut_int) C Calculate the angle-dependent terms of energy & contributions to derivatives. call sc_angular C Calculate whole angle-dependent part of epsilon and contributions @@ -480,7 +534,7 @@ C to its derivatives eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij*(1.0d0-sss) + evdw=evdw+evdwij*(1.0d0-sss)*sss1 if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa @@ -496,13 +550,15 @@ C Calculate gradient components. fac=-expon*(e1+evdwij) sigder=fac/sigsq fac=rrij*fac + & +evdwij*(sss1*sssgrad/sigmaii(itypi,itypj) + & +(1.0d0-sss)*sssgrad1)*rij C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac 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) + call sc_grad_scale((1.0d0-sss)*sss1) endif enddo ! j enddo ! iint @@ -527,7 +583,13 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include "COMMON.SPLITELE" + integer icall common /srutu/ icall + double precision evdw + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi + double precision sscale,sscagrad c double precision rrsave(maxdim) logical lprn evdw=0.0D0 @@ -540,9 +602,9 @@ c else c endif ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -557,7 +619,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -578,7 +640,7 @@ c dscj_inv=dsc_inv(itypj) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) + sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) if (sss.gt.0.0d0) then @@ -609,6 +671,7 @@ C Calculate gradient components. fac=-expon*(e1+evdwij) sigder=fac/sigsq fac=rrij*fac + & +evdwij*sssgrad/sigmaii(itypi,itypj)*rrij C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac @@ -629,7 +692,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -641,8 +704,17 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include "COMMON.SPLITELE" logical lprn integer xshift,yshift,zshift + double precision evdw + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi + double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, + & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + double precision dist,sscale,sscagrad,sscagradlip,sscalelip + double precision subchap,sss1,sssgrad1 evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -651,9 +723,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.false. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -676,7 +748,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -696,81 +768,81 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) xj=c(1,nres+j) yj=c(2,nres+j) zj=c(3,nres+j) - xj=mod(xj,boxxsize) - if (xj.lt.0) xj=xj+boxxsize - yj=mod(yj,boxysize) - if (yj.lt.0) yj=yj+boxysize - zj=mod(zj,boxzsize) - if (zj.lt.0) zj=zj+boxzsize - if ((zj.gt.bordlipbot) - &.and.(zj.lt.bordliptop)) then + 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 + if (zj.lt.buflipbot) then C what fraction I am in - fracinbuf=1.0d0- - & ((positi-bordlipbot)/lipbufthick) + fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore - sslipj=sscalelip(fracinbuf) - ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) - sslipj=sscalelip(fracinbuf) - ssgradlipj=sscagradlip(fracinbuf)/lipbufthick - else - sslipj=1.0d0 - ssgradlipj=0.0 - endif - else - sslipj=0.0d0 - ssgradlipj=0.0 - endif - aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 - bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 - - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 - xj=xj_safe+xshift*boxxsize - yj=yj_safe+yshift*boxysize - zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - subchap=1 - endif - enddo - enddo - enddo - if (subchap.eq.1) then - xj=xj_temp-xi - yj=yj_temp-yi - zj=zj_temp-zi - else - xj=xj_safe-xi - yj=yj_safe-yi - zj=zj_safe-zi - endif - dxj=dc_norm(1,nres+j) - dyj=dc_norm(2,nres+j) - dzj=dc_norm(3,nres+j) - rrij=1.0D0/(xj*xj+yj*yj+zj*zj) - rij=dsqrt(rrij) - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) - sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj)) - if (sss.lt.1.0d0) then - + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + else if (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if (dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif + dxj=dc_norm(1,nres+j) + dyj=dc_norm(2,nres+j) + dzj=dc_norm(3,nres+j) + rrij=1.0D0/(xj*xj+yj*yj+zj*zj) + rij=dsqrt(rrij) + sss1=sscale(1.0d0/rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + 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. + sssgrad= + & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) + sssgrad1=sscagrad(1.0d0/rij,r_cut_int) call sc_angular sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) @@ -797,7 +869,7 @@ c--------------------------------------------------------------- c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij*(1.0d0-sss) + evdw=evdw+evdwij*(1.0d0-sss)*sss1 if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa @@ -809,15 +881,16 @@ 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,'(a6,2i5,4f10.5)') + & 'evdw',i,j,rij,sss,sss1,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac - fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij + fac=fac+evdwij*(-sss1*sssgrad/sigmaii(itypi,itypj) + & +(1.0d0-sss)*sssgrad1)*rij c fac=0.0d0 C Calculate the radial part of the gradient gg(1)=xj*fac @@ -826,7 +899,7 @@ C Calculate the radial part of the gradient gg_lipi(3)=ssgradlipi*evdwij gg_lipj(3)=ssgradlipj*evdwij C Calculate angular part of the gradient. - call sc_grad_scale(1.0d0-sss) + call sc_grad_scale((1.0d0-sss)*sss1) endif enddo ! j enddo ! iint @@ -841,7 +914,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -853,8 +926,17 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include "COMMON.SPLITELE" logical lprn integer xshift,yshift,zshift + double precision evdw + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi + double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, + & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + double precision dist,sscale,sscagrad,sscagradlip,sscalelip + double precision subchap evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -863,9 +945,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.false. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -888,7 +970,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -908,76 +990,74 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) xj=c(1,nres+j) yj=c(2,nres+j) zj=c(3,nres+j) - xj=mod(xj,boxxsize) - if (xj.lt.0) xj=xj+boxxsize - yj=mod(yj,boxysize) - if (yj.lt.0) yj=yj+boxysize - zj=mod(zj,boxzsize) - if (zj.lt.0) zj=zj+boxzsize - if ((zj.gt.bordlipbot) - &.and.(zj.lt.bordliptop)) then + 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 + if (zj.lt.buflipbot) then C what fraction I am in - fracinbuf=1.0d0- - & ((positi-bordlipbot)/lipbufthick) + fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore - sslipj=sscalelip(fracinbuf) - ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) - sslipj=sscalelip(fracinbuf) - ssgradlipj=sscagradlip(fracinbuf)/lipbufthick - else - sslipj=1.0d0 - ssgradlipj=0.0 - endif - else - sslipj=0.0d0 - ssgradlipj=0.0 - endif - aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 - bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 - xj=xj_safe+xshift*boxxsize - yj=yj_safe+yshift*boxysize - zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - subchap=1 - endif - enddo - enddo - enddo - if (subchap.eq.1) then - xj=xj_temp-xi - yj=yj_temp-yi - zj=zj_temp-zi - else - xj=xj_safe-xi - yj=yj_safe-yi - zj=zj_safe-zi - endif + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) - sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj)) + 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 C Calculate angle-dependent terms of energy and contributions to their @@ -1028,7 +1108,7 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac - fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij + fac=fac+evdwij*sssgrad/sigmaii(itypi,itypj)*rij c fac=0.0d0 C Calculate the radial part of the gradient gg(1)=xj*fac @@ -1063,8 +1143,18 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include "COMMON.SPLITELE" + integer icall common /srutu/ icall logical lprn + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, + & xi,yi,zi,fac_augm,e_augm + double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, + & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + double precision dist,sscale,sscagrad,sscagradlip,sscalelip + double precision sss1,sssgrad1 evdw=0.0D0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 @@ -1072,9 +1162,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1089,7 +1179,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -1113,9 +1203,13 @@ c dscj_inv=dsc_inv(itypj) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) + sss1=sscale(1.0d0/rij,r_cut_int) + if (sss1.eq.0.0d0) cycle if (sss.lt.1.0d0) then + sssgrad= + & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) + sssgrad1=sscagrad(1.0d0/rij,r_cut_int) C Calculate angle-dependent terms of energy and contributions to their C derivatives. @@ -1140,7 +1234,7 @@ c--------------------------------------------------------------- fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+(evdwij+e_augm)*(1.0d0-sss) + evdw=evdw+(evdwij+e_augm)*sss1*(1.0d0-sss) if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa @@ -1157,12 +1251,15 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac-2*expon*rrij*e_augm + fac=fac+(evdwij+e_augm)* + & (-sss1*sssgrad/sigmaii(itypi,itypj) + & +(1.0d0-sss)*sssgrad1)*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac C Calculate angular part of the gradient. - call sc_grad_scale(1.0d0-sss) + call sc_grad_scale((1.0d0-sss)*sss1) endif enddo ! j enddo ! iint @@ -1174,7 +1271,7 @@ C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne-Vorobjev potential of interaction. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -1185,8 +1282,18 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include "COMMON.SPLITELE" + integer icall common /srutu/ icall logical lprn + integer itypi,itypj,itypi1,iint,ind + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, + & xi,yi,zi,fac_augm,e_augm + double precision evdw + double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, + & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + double precision dist,sscale,sscagrad,sscagradlip,sscalelip evdw=0.0D0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 @@ -1194,9 +1301,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1211,7 +1318,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -1235,7 +1342,7 @@ c dscj_inv=dsc_inv(itypj) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) + sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) if (sss.gt.0.0d0) then @@ -1278,7 +1385,8 @@ C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder - fac=rij*fac-2*expon*rrij*e_augm + fac=rij*fac-2*expon*rrij*e_augm+ + & (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac @@ -1298,6 +1406,7 @@ C---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.CALC' include 'COMMON.IOUNITS' + include "COMMON.SPLITELE" double precision dcosom1(3),dcosom2(3) double precision scalfac eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 @@ -1365,12 +1474,17 @@ C include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' +#ifdef FOURBODY include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' +#endif + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' include 'COMMON.SHIELD' + include "COMMON.SPLITELE" dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -1439,9 +1553,11 @@ cd enddo eello_turn3=0.0d0 eello_turn4=0.0d0 ind=0 +#ifdef FOURBODY do i=1,nres num_cont_hb(i)=0 enddo +#endif cd print '(a)','Enter EELEC' cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e do i=1,nres @@ -1478,7 +1594,9 @@ C & .or. itype(i+4).eq.ntyp1 num_conti=0 call eelecij_scale(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo do i=iturn4_start,iturn4_end if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 @@ -1502,11 +1620,15 @@ C & .or. itype(i-1).eq.ntyp1 if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize +#ifdef FOURBODY num_conti=num_cont_hb(i) +#endif call eelecij_scale(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo ! i c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 @@ -1531,8 +1653,10 @@ C & .or. itype(i-1).eq.ntyp1 if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize -c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) +c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend +#ifdef FOURBODY num_conti=num_cont_hb(i) +#endif do j=ielstart(i),ielend(i) if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 C & .or.itype(j+2).eq.ntyp1 @@ -1540,7 +1664,9 @@ C & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij_scale(i,j,ees,evdw1,eel_loc) enddo ! j +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo ! i c write (iout,*) "Number of loop steps in EELEC:",ind cd do i=1,nres @@ -1554,7 +1680,7 @@ cd print *,"Processor",fg_rank," t_eelecij",t_eelecij end C------------------------------------------------------------------------------- subroutine eelecij_scale(i,j,ees,evdw1,eel_loc) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" @@ -1567,21 +1693,48 @@ C------------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' +#ifdef FOURBODY include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' +#endif + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' include 'COMMON.SHIELD' + include "COMMON.SPLITELE" integer xshift,yshift,zshift - dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), + double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4), & gmuij2(4),gmuji2(4) + integer j1,j2,num_conti 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 + integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj + integer ilist,iresshield + double precision rlocshield + double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp + double precision ees,evdw1,eel_loc,aaa,bbb,ael3i + double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj, + & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4, + & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa, + & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der, + & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij, + & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp, + & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp, + & ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji, + & dxi,dyi,dzi,a22,a23,a32,a33 + double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe, + & zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi, + & dz_normi,aux + double precision sss1,sssgrad1 + double precision sscale,sscagrad + double precision scalar + c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -1596,29 +1749,29 @@ C 13-go grudnia roku pamietnego... c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j C print *,"WCHODZE2" - ind=ind+1 - iteli=itel(i) - itelj=itel(j) - if (j.eq.i+2 .and. itelj.eq.2) iteli=2 - aaa=app(iteli,itelj) - bbb=bpp(iteli,itelj) - ael6i=ael6(iteli,itelj) - ael3i=ael3(iteli,itelj) - dxj=dc(1,j) - dyj=dc(2,j) - dzj=dc(3,j) - 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 - 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 + ind=ind+1 + iteli=itel(i) + itelj=itel(j) + if (j.eq.i+2 .and. itelj.eq.2) iteli=2 + aaa=app(iteli,itelj) + bbb=bpp(iteli,itelj) + ael6i=ael6(iteli,itelj) + ael3i=ael3(iteli,itelj) + dxj=dc(1,j) + dyj=dc(2,j) + dzj=dc(3,j) + 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 + 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 @@ -1638,89 +1791,102 @@ C print *,"WCHODZE2" 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 + 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 - rrmij=1.0D0/rij - rij=dsqrt(rij) - rmij=1.0D0/rij + rij=xj*xj+yj*yj+zj*zj + rrmij=1.0D0/rij + rij=dsqrt(rij) + rmij=1.0D0/rij c For extracting the short-range part of Evdwpp - sss=sscale(rij/rpp(iteli,itelj)) - sssgrad=sscagrad(rij/rpp(iteli,itelj)) - r3ij=rrmij*rmij - r6ij=r3ij*r3ij - cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj - cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij - cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij - fac=cosa-3.0D0*cosb*cosg - ev1=aaa*r6ij*r6ij + sss1=sscale(rij,r_cut_int) + if (sss1.eq.0.0d0) return + sss=sscale(rij/rpp(iteli,itelj),r_cut_respa) + sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa) + sssgrad1=sscagrad(rij,r_cut_int) + r3ij=rrmij*rmij + r6ij=r3ij*r3ij + cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj + cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij + cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij + fac=cosa-3.0D0*cosb*cosg + ev1=aaa*r6ij*r6ij c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions - if (j.eq.i+2) ev1=scal_el*ev1 - ev2=bbb*r6ij - fac3=ael6i*r6ij - fac4=ael3i*r3ij - evdwij=ev1+ev2 - if (shield_mode.eq.0) then - fac_shield(i)=1.0 - fac_shield(j)=1.0 - endif - el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) - el2=fac4*fac - el1=el1*fac_shield(i)**2*fac_shield(j)**2 - el2=el2*fac_shield(i)**2*fac_shield(j)**2 - eesij=el1+el2 + if (j.eq.i+2) ev1=scal_el*ev1 + ev2=bbb*r6ij + fac3=ael6i*r6ij + fac4=ael3i*r3ij + evdwij=ev1+ev2 + if (shield_mode.eq.0) then + fac_shield(i)=1.0 + fac_shield(j)=1.0 + endif + el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) + el2=fac4*fac + el1=el1*fac_shield(i)**2*fac_shield(j)**2 + el2=el2*fac_shield(i)**2*fac_shield(j)**2 + 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 - evdw1=evdw1+evdwij*(1.0d0-sss) + ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) + ees=ees+eesij*sss1 + evdw1=evdw1+evdwij*(1.0d0-sss)*sss1 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,f7.3)') 'evdw1',i,j,evdwij,sss - write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij - endif + 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 + endif C C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE - facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss) - facel=-3*rrmij*(el1+eesij) - fac1=fac - erij(1)=xj*rmij - erij(2)=yj*rmij - erij(3)=zj*rmij + facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1 +c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + facel=-3*rrmij*(el1+eesij)*sss1 +c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + fac1=fac + erij(1)=xj*rmij + erij(2)=yj*rmij + erij(3)=zj*rmij * * Radial derivatives. First process both termini of the fragment (i,j) * - ggg(1)=facel*xj - ggg(2)=facel*yj - ggg(3)=facel*zj - if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + aux=facel+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj)) + & *eesij*rmij +c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=aux*zj +c ggg(1)=facel*xj +c ggg(2)=facel*yj +c ggg(3)=facel*zj + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j - do ilist=1,ishield_list(i) - iresshield=shield_list(ilist,i) - do k=1,3 - rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i) - & *2.0 + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1 + & /fac_shield(i)*2.0*sss1 gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0 + & +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0 + & *sss1 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) @@ -1737,32 +1903,32 @@ C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) C C enddo C endif - enddo enddo - do ilist=1,ishield_list(j) - iresshield=shield_list(ilist,j) - do k=1,3 + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) - & *2.0 + & *2.0*sss1 gshieldx(k,iresshield)=gshieldx(k,iresshield)+ - & rlocshield - & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 + & rlocshield + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield - enddo enddo + enddo - do k=1,3 - gshieldc(k,i)=gshieldc(k,i)+ - & grad_shield(k,i)*eesij/fac_shield(i)*2.0 - gshieldc(k,j)=gshieldc(k,j)+ - & grad_shield(k,j)*eesij/fac_shield(j)*2.0 - gshieldc(k,i-1)=gshieldc(k,i-1)+ - & grad_shield(k,i)*eesij/fac_shield(i)*2.0 - gshieldc(k,j-1)=gshieldc(k,j-1)+ - & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + do k=1,3 + gshieldc(k,i)=gshieldc(k,i)+ + & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1 + gshieldc(k,j)=gshieldc(k,j)+ + & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1 + gshieldc(k,i-1)=gshieldc(k,i-1)+ + & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1 + gshieldc(k,j-1)=gshieldc(k,j-1)+ + & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1 - enddo - endif + enddo + endif c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -1770,10 +1936,12 @@ c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end - do k=1,3 - gelc_long(k,j)=gelc_long(k,j)+ggg(k) - gelc_long(k,i)=gelc_long(k,i)-ggg(k) - enddo + do k=1,3 + gelc_long(k,j)=gelc_long(k,j)+ggg(k) + gelc_long(k,i)=gelc_long(k,i)-ggg(k) + enddo +c gelc_long(3,i)=gelc_long(3,i)+ +c ssgradlipi*eesij/2.0d0*lipscale**2*sss1 * * Loop over residues i+1 thru j-1. * @@ -1782,19 +1950,22 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - 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=facvdw+ + & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij +c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(1)=facvdw*xj + ggg(2)=facvdw*yj + ggg(3)=facvdw*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf c gvdwpp(k,j)=gvdwpp(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end - do k=1,3 - gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) - gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) - enddo + do k=1,3 + gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) + gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) + enddo * * Loop over residues i+1 thru j-1. * @@ -1804,29 +1975,40 @@ cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) cgrad enddo cgrad enddo #else - facvdw=ev1+evdwij*(1.0d0-sss) - facel=el1+eesij - fac1=fac - fac=-3*rrmij*(facvdw+facvdw+facel) - erij(1)=xj*rmij - erij(2)=yj*rmij - erij(3)=zj*rmij + facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1 +c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + facel=-3*rrmij*(el1+eesij)*sss1 +c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + +c facvdw=ev1+evdwij*(1.0d0-sss)*sss1 +c facel=el1+eesij + fac1=fac + fac=-3*rrmij*(facvdw+facvdw+facel) + erij(1)=xj*rmij + erij(2)=yj*rmij + erij(3)=zj*rmij * * Radial derivatives. First process both termini of the fragment (i,j) * - ggg(1)=fac*xj - ggg(2)=fac*yj - ggg(3)=fac*zj + aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj)) + & *eesij*rmij +c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=axu*zj +c ggg(1)=fac*xj +c ggg(2)=fac*yj +c ggg(3)=fac*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end - do k=1,3 - gelc_long(k,j)=gelc(k,j)+ggg(k) - gelc_long(k,i)=gelc(k,i)-ggg(k) - enddo + do k=1,3 + gelc_long(k,j)=gelc(k,j)+ggg(k) + gelc_long(k,i)=gelc(k,i)-ggg(k) + enddo * * Loop over residues i+1 thru j-1. * @@ -1839,33 +2021,36 @@ c 9/28/08 AL Gradient compotents will be summed only at the end C ggg(1)=facvdw*xj C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj - 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) - do k=1,3 - gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) - gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) - enddo + facvdw=facvdw + & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij + ggg(1)=facvdw*xj + ggg(2)=facvdw*yj + ggg(3)=facvdw*zj + do k=1,3 + gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) + gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) + enddo #endif * * Angular part * - ecosa=2.0D0*fac3*fac1+fac4 - fac4=-3.0D0*fac4 - fac3=-6.0D0*fac3 - ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4) - ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4) - do k=1,3 - dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb) - dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg) - enddo + ecosa=2.0D0*fac3*fac1+fac4 + fac4=-3.0D0*fac4 + fac3=-6.0D0*fac3 + ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4) + ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4) + do k=1,3 + dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb) + dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg) + enddo 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))* - & fac_shield(i)**2*fac_shield(j)**2 + do k=1,3 + ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1 + & *fac_shield(i)**2*fac_shield(j)**2 +c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) - enddo + enddo c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf @@ -1880,22 +2065,24 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - do k=1,3 - 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)) - & *fac_shield(i)**2*fac_shield(j)**2 + do k=1,3 + 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 +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)) - & *fac_shield(i)**2*fac_shield(j)**2 - gelc_long(k,j)=gelc_long(k,j)+ggg(k) - gelc_long(k,i)=gelc_long(k,i)-ggg(k) - enddo - IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 - & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 - & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN + 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 +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) + enddo + IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 + & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 + & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN C C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction C energy of a peptide unit is assumed in the form of a second-order @@ -1903,44 +2090,44 @@ C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al. C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms C are computed for EVERY pair of non-contiguous peptide groups. C - if (j.lt.nres-1) then - j1=j+1 - j2=j-1 - else - j1=j-1 - j2=j-2 - endif - kkk=0 - do k=1,2 - do l=1,2 - kkk=kkk+1 - muij(kkk)=mu(k,i)*mu(l,j) + if (j.lt.nres-1) then + j1=j+1 + j2=j-1 + else + j1=j-1 + j2=j-2 + endif + kkk=0 + do k=1,2 + do l=1,2 + kkk=kkk+1 + muij(kkk)=mu(k,i)*mu(l,j) #ifdef NEWCORR - gmuij1(kkk)=gtb1(k,i+1)*mu(l,j) + gmuij1(kkk)=gtb1(k,i+1)*mu(l,j) c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j) - gmuij2(kkk)=gUb2(k,i)*mu(l,j) - gmuji1(kkk)=mu(k,i)*gtb1(l,j+1) + gmuij2(kkk)=gUb2(k,i)*mu(l,j) + gmuji1(kkk)=mu(k,i)*gtb1(l,j+1) c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i) - gmuji2(kkk)=mu(k,i)*gUb2(l,j) + gmuji2(kkk)=mu(k,i)*gUb2(l,j) #endif - enddo - enddo + enddo + enddo cd write (iout,*) 'EELEC: i',i,' j',j cd write (iout,*) 'j',j,' j1',j1,' j2',j2 cd write(iout,*) 'muij',muij - ury=scalar(uy(1,i),erij) - urz=scalar(uz(1,i),erij) - vry=scalar(uy(1,j),erij) - vrz=scalar(uz(1,j),erij) - a22=scalar(uy(1,i),uy(1,j))-3*ury*vry - a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz - a32=scalar(uz(1,i),uy(1,j))-3*urz*vry - a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz - fac=dsqrt(-ael6i)*r3ij - a22=a22*fac - a23=a23*fac - a32=a32*fac - a33=a33*fac + ury=scalar(uy(1,i),erij) + urz=scalar(uz(1,i),erij) + vry=scalar(uy(1,j),erij) + vrz=scalar(uz(1,j),erij) + a22=scalar(uy(1,i),uy(1,j))-3*ury*vry + a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz + a32=scalar(uz(1,i),uy(1,j))-3*urz*vry + a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz + fac=dsqrt(-ael6i)*r3ij + a22=a22*fac + a23=a23*fac + a32=a32*fac + a33=a33*fac cd write (iout,'(4i5,4f10.5)') cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij @@ -1953,101 +2140,113 @@ cd write (iout,'(4f10.5)') ury,urz,vry,vrz cd write (iout,'(9f10.5/)') cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij C Derivatives of the elements of A in virtual-bond vectors - call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1)) - do k=1,3 - uryg(k,1)=scalar(erder(1,k),uy(1,i)) - uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1)) - uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1)) - urzg(k,1)=scalar(erder(1,k),uz(1,i)) - urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1)) - urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1)) - vryg(k,1)=scalar(erder(1,k),uy(1,j)) - vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1)) - vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1)) - vrzg(k,1)=scalar(erder(1,k),uz(1,j)) - vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1)) - vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1)) - enddo + call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1)) + do k=1,3 + uryg(k,1)=scalar(erder(1,k),uy(1,i)) + uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1)) + uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1)) + urzg(k,1)=scalar(erder(1,k),uz(1,i)) + urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1)) + urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1)) + vryg(k,1)=scalar(erder(1,k),uy(1,j)) + vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1)) + vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1)) + vrzg(k,1)=scalar(erder(1,k),uz(1,j)) + vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1)) + vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1)) + enddo C Compute radial contributions to the gradient - facr=-3.0d0*rrmij - a22der=a22*facr - a23der=a23*facr - a32der=a32*facr - a33der=a33*facr - agg(1,1)=a22der*xj - agg(2,1)=a22der*yj - agg(3,1)=a22der*zj - agg(1,2)=a23der*xj - agg(2,2)=a23der*yj - agg(3,2)=a23der*zj - agg(1,3)=a32der*xj - agg(2,3)=a32der*yj - agg(3,3)=a32der*zj - agg(1,4)=a33der*xj - agg(2,4)=a33der*yj - agg(3,4)=a33der*zj + facr=-3.0d0*rrmij + a22der=a22*facr + a23der=a23*facr + a32der=a32*facr + a33der=a33*facr + agg(1,1)=a22der*xj + agg(2,1)=a22der*yj + agg(3,1)=a22der*zj + agg(1,2)=a23der*xj + agg(2,2)=a23der*yj + agg(3,2)=a23der*zj + agg(1,3)=a32der*xj + agg(2,3)=a32der*yj + agg(3,3)=a32der*zj + agg(1,4)=a33der*xj + agg(2,4)=a33der*yj + agg(3,4)=a33der*zj C Add the contributions coming from er - fac3=-3.0d0*fac - do k=1,3 - agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury) - agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury) - agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz) - agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz) - enddo - do k=1,3 + fac3=-3.0d0*fac + do k=1,3 + agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury) + agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury) + agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz) + agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz) + enddo + do k=1,3 C Derivatives in DC(i) cgrad ghalf1=0.5d0*agg(k,1) cgrad ghalf2=0.5d0*agg(k,2) cgrad ghalf3=0.5d0*agg(k,3) cgrad ghalf4=0.5d0*agg(k,4) - aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j)) - & -3.0d0*uryg(k,2)*vry)!+ghalf1 - aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j)) - & -3.0d0*uryg(k,2)*vrz)!+ghalf2 - aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j)) - & -3.0d0*urzg(k,2)*vry)!+ghalf3 - aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j)) - & -3.0d0*urzg(k,2)*vrz)!+ghalf4 + aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j)) + & -3.0d0*uryg(k,2)*vry)!+ghalf1 + aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j)) + & -3.0d0*uryg(k,2)*vrz)!+ghalf2 + aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j)) + & -3.0d0*urzg(k,2)*vry)!+ghalf3 + aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j)) + & -3.0d0*urzg(k,2)*vrz)!+ghalf4 C Derivatives in DC(i+1) - aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j)) - & -3.0d0*uryg(k,3)*vry)!+agg(k,1) - aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j)) - & -3.0d0*uryg(k,3)*vrz)!+agg(k,2) - aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j)) - & -3.0d0*urzg(k,3)*vry)!+agg(k,3) - aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j)) - & -3.0d0*urzg(k,3)*vrz)!+agg(k,4) + aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j)) + & -3.0d0*uryg(k,3)*vry)!+agg(k,1) + aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j)) + & -3.0d0*uryg(k,3)*vrz)!+agg(k,2) + aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j)) + & -3.0d0*urzg(k,3)*vry)!+agg(k,3) + aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j)) + & -3.0d0*urzg(k,3)*vrz)!+agg(k,4) C Derivatives in DC(j) - aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i)) - & -3.0d0*vryg(k,2)*ury)!+ghalf1 - aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i)) - & -3.0d0*vrzg(k,2)*ury)!+ghalf2 - aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i)) - & -3.0d0*vryg(k,2)*urz)!+ghalf3 - aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i)) - & -3.0d0*vrzg(k,2)*urz)!+ghalf4 + aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i)) + & -3.0d0*vryg(k,2)*ury)!+ghalf1 + aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i)) + & -3.0d0*vrzg(k,2)*ury)!+ghalf2 + aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i)) + & -3.0d0*vryg(k,2)*urz)!+ghalf3 + aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i)) + & -3.0d0*vrzg(k,2)*urz)!+ghalf4 C Derivatives in DC(j+1) or DC(nres-1) - aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i)) - & -3.0d0*vryg(k,3)*ury) - aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i)) - & -3.0d0*vrzg(k,3)*ury) - aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i)) - & -3.0d0*vryg(k,3)*urz) - aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i)) - & -3.0d0*vrzg(k,3)*urz) + aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i)) + & -3.0d0*vryg(k,3)*ury) + aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i)) + & -3.0d0*vrzg(k,3)*ury) + aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i)) + & -3.0d0*vryg(k,3)*urz) + aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i)) + & -3.0d0*vrzg(k,3)*urz) cgrad if (j.eq.nres-1 .and. i.lt.j-2) then cgrad do l=1,4 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l) cgrad enddo cgrad endif + enddo + acipa(1,1)=a22 + acipa(1,2)=a23 + acipa(2,1)=a32 + acipa(2,2)=a33 + a22=-a22 + a23=-a23 + do l=1,2 + do k=1,3 + agg(k,l)=-agg(k,l) + aggi(k,l)=-aggi(k,l) + aggi1(k,l)=-aggi1(k,l) + aggj(k,l)=-aggj(k,l) + aggj1(k,l)=-aggj1(k,l) enddo - acipa(1,1)=a22 - acipa(1,2)=a23 - acipa(2,1)=a32 - acipa(2,2)=a33 + enddo + if (j.lt.nres-1) then a22=-a22 - a23=-a23 - do l=1,2 + a32=-a32 + do l=1,3,2 do k=1,3 agg(k,l)=-agg(k,l) aggi(k,l)=-aggi(k,l) @@ -2056,56 +2255,44 @@ cgrad endif aggj1(k,l)=-aggj1(k,l) enddo enddo - if (j.lt.nres-1) then - a22=-a22 - a32=-a32 - do l=1,3,2 - do k=1,3 - agg(k,l)=-agg(k,l) - aggi(k,l)=-aggi(k,l) - aggi1(k,l)=-aggi1(k,l) - aggj(k,l)=-aggj(k,l) - aggj1(k,l)=-aggj1(k,l) - enddo + else + a22=-a22 + a23=-a23 + a32=-a32 + a33=-a33 + do l=1,4 + do k=1,3 + agg(k,l)=-agg(k,l) + aggi(k,l)=-aggi(k,l) + aggi1(k,l)=-aggi1(k,l) + aggj(k,l)=-aggj(k,l) + aggj1(k,l)=-aggj1(k,l) enddo - else - a22=-a22 - a23=-a23 - a32=-a32 - a33=-a33 - do l=1,4 - do k=1,3 - agg(k,l)=-agg(k,l) - aggi(k,l)=-aggi(k,l) - aggi1(k,l)=-aggi1(k,l) - aggj(k,l)=-aggj(k,l) - aggj1(k,l)=-aggj1(k,l) - enddo - enddo - endif - ENDIF ! WCORR - IF (wel_loc.gt.0.0d0) THEN + enddo + endif + ENDIF ! WCORR + IF (wel_loc.gt.0.0d0) THEN 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 + 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 - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'eelloc',i,j,eel_loc_ij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') + & 'eelloc',i,j,eel_loc_ij - if (shield_mode.eq.0) then - fac_shield(i)=1.0 - fac_shield(j)=1.0 -C else -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) - eel_loc=eel_loc+eel_loc_ij + if (shield_mode.eq.0) then + fac_shield(i)=1.0 + fac_shield(j)=1.0 +C else +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) + eel_loc=eel_loc+eel_loc_ij*sss1 - if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j @@ -2113,11 +2300,11 @@ C print *,i,j iresshield=shield_list(ilist,i) do k=1,3 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij - & /fac_shield(i) + & /fac_shield(i)*sss1 C & *2.0 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i) + & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)*sss1 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) & +rlocshield enddo @@ -2126,11 +2313,11 @@ C & *2.0 iresshield=shield_list(ilist,j) do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij - & /fac_shield(j) + & /fac_shield(j)*sss1 C & *2.0 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j) + & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)*sss1 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) & +rlocshield @@ -2139,41 +2326,41 @@ C & *2.0 do k=1,3 gshieldc_ll(k,i)=gshieldc_ll(k,i)+ - & grad_shield(k,i)*eel_loc_ij/fac_shield(i) + & grad_shield(k,i)*eel_loc_ij/fac_shield(i)*sss1 gshieldc_ll(k,j)=gshieldc_ll(k,j)+ - & grad_shield(k,j)*eel_loc_ij/fac_shield(j) + & grad_shield(k,j)*eel_loc_ij/fac_shield(j)*sss1 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+ - & grad_shield(k,i)*eel_loc_ij/fac_shield(i) + & grad_shield(k,i)*eel_loc_ij/fac_shield(i)*sss1 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+ - & grad_shield(k,j)*eel_loc_ij/fac_shield(j) - enddo - endif + & grad_shield(k,j)*eel_loc_ij/fac_shield(j)*sss1 + enddo + endif #ifdef NEWCORR - geel_loc_ij=(a22*gmuij1(1) + geel_loc_ij=(a22*gmuij1(1) & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss1 c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) - gloc(nphi+i,icg)=gloc(nphi+i,icg)+ + gloc(nphi+i,icg)=gloc(nphi+i,icg)+ & geel_loc_ij*wel_loc c write(iout,*) "derivative over thatai-1" c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), c & a33*gmuij2(4) - geel_loc_ij= + geel_loc_ij= & a22*gmuij2(1) & +a23*gmuij2(2) & +a32*gmuij2(3) & +a33*gmuij2(4) - gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ + gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss1 c Derivative over j residue - geel_loc_ji=a22*gmuji1(1) + geel_loc_ji=a22*gmuji1(1) & +a23*gmuji1(2) & +a32*gmuji1(3) & +a33*gmuji1(4) @@ -2183,9 +2370,9 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss1 - geel_loc_ji= + geel_loc_ji= & +a22*gmuji2(1) & +a23*gmuji2(2) & +a32*gmuji2(3) @@ -2193,147 +2380,167 @@ c & a33*gmuji1(4) c write(iout,*) "derivative over thataj-1" 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) + gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ + & geel_loc_ji*wel_loc + & *fac_shield(i)*fac_shield(j)*sss1 #endif -cC Partial derivatives in virtual-bond dihedral angles gamma - if (i.gt.1) +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) + & (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 +c & *fac_shield(i)*fac_shield(j) +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) + + 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 +c & *fac_shield(i)*fac_shield(j) +c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) - do l=1,3 - 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) + do l=1,3 + 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 +c & *fac_shield(i)*fac_shield(j) +c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) - gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) - gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) -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(l,j)=gel_loc_long(l,j)+ggg(l) + gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) +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 cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) cgrad enddo cgrad enddo C Remaining derivatives of eello - 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) +c gel_loc_long(3,j)=gel_loc_long(3,j)+ & +c ssgradlipj*eel_loc_ij/2.0d0*lipscale/ & +c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut +c +c gel_loc_long(3,i)=gel_loc_long(3,i)+ & +c ssgradlipi*eel_loc_ij/2.0d0*lipscale/ & +c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut - 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) + 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 +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) + 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 +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) + 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 +c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) - enddo - ENDIF + 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 +c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + + enddo + ENDIF +#ifdef FOURBODY C Change 12/26/95 to calculate four-body contributions to H-bonding energy c if (j.gt.i+1 .and. num_conti.le.maxconts) then - if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0 - & .and. num_conti.le.maxconts) then + if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0 + & .and. num_conti.le.maxconts) then c write (iout,*) i,j," entered corr" C C Calculate the contact function. The ith column of the array JCONT will C contain the numbers of atoms that make contacts with the atom I (of numbers C greater than I). The arrays FACONT and GACONT will contain the values of C the contact function and its derivative. -c r0ij=1.02D0*rpp(iteli,itelj) -c r0ij=1.11D0*rpp(iteli,itelj) - r0ij=2.20D0*rpp(iteli,itelj) -c r0ij=1.55D0*rpp(iteli,itelj) - call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont) - if (fcont.gt.0.0D0) then - num_conti=num_conti+1 - if (num_conti.gt.maxconts) then - write (iout,*) 'WARNING - max. # of contacts exceeded;', - & ' will skip next contacts for this conf.' - else - jcont_hb(num_conti,i)=j -cd write (iout,*) "i",i," j",j," num_conti",num_conti, -cd & " jcont_hb",jcont_hb(num_conti,i) - IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. - & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN +c r0ij=1.02D0*rpp(iteli,itelj) +c r0ij=1.11D0*rpp(iteli,itelj) + r0ij=2.20D0*rpp(iteli,itelj) +c r0ij=1.55D0*rpp(iteli,itelj) + call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont) + if (fcont.gt.0.0D0) then + num_conti=num_conti+1 + if (num_conti.gt.maxconts) then + write (iout,*) 'WARNING - max. # of contacts exceeded;', + & ' will skip next contacts for this conf.' + else + jcont_hb(num_conti,i)=j +cd write (iout,*) "i",i," j",j," num_conti",num_conti, +cd " jcont_hb",jcont_hb(num_conti,i) + IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. + & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el C terms. - d_cont(num_conti,i)=rij + d_cont(num_conti,i)=rij cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij C --- Electrostatic-interaction matrix --- - a_chuj(1,1,num_conti,i)=a22 - a_chuj(1,2,num_conti,i)=a23 - a_chuj(2,1,num_conti,i)=a32 - a_chuj(2,2,num_conti,i)=a33 + a_chuj(1,1,num_conti,i)=a22 + a_chuj(1,2,num_conti,i)=a23 + a_chuj(2,1,num_conti,i)=a32 + a_chuj(2,2,num_conti,i)=a33 C --- Gradient of rij - do kkk=1,3 - grij_hb_cont(kkk,num_conti,i)=erij(kkk) - enddo - kkll=0 - do k=1,2 - do l=1,2 - kkll=kkll+1 - do m=1,3 - a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll) - a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll) - a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll) - a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll) - a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll) - enddo + do kkk=1,3 + grij_hb_cont(kkk,num_conti,i)=erij(kkk) + enddo + kkll=0 + do k=1,2 + do l=1,2 + kkll=kkll+1 + do m=1,3 + a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll) + a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll) + a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll) + a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll) + a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll) enddo enddo - ENDIF - IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN + enddo + ENDIF + IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN C Calculate contact energies - cosa4=4.0D0*cosa - wij=cosa-3.0D0*cosb*cosg - cosbg1=cosb+cosg - cosbg2=cosb-cosg -c fac3=dsqrt(-ael6i)/r0ij**3 - fac3=dsqrt(-ael6i)*r3ij -c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1) - ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1 - if (ees0tmp.gt.0) then - ees0pij=dsqrt(ees0tmp) - else - ees0pij=0 - endif -c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) - ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2 - if (ees0tmp.gt.0) then - ees0mij=dsqrt(ees0tmp) - else - ees0mij=0 - endif -c ees0mij=0.0D0 - if (shield_mode.eq.0) then + cosa4=4.0D0*cosa + wij=cosa-3.0D0*cosb*cosg + cosbg1=cosb+cosg + cosbg2=cosb-cosg +c fac3=dsqrt(-ael6i)/r0ij**3 + fac3=dsqrt(-ael6i)*r3ij +c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1) + ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1 + if (ees0tmp.gt.0) then + ees0pij=dsqrt(ees0tmp) + else + ees0pij=0 + endif +c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) + ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2 + if (ees0tmp.gt.0) then + ees0mij=dsqrt(ees0tmp) + else + ees0mij=0 + endif +c ees0mij=0.0D0 + if (shield_mode.eq.0) then fac_shield(i)=1.0d0 fac_shield(j)=1.0d0 - else + else ees0plist(num_conti,i)=j C fac_shield(i)=0.4d0 C fac_shield(j)=0.6d0 - endif - ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) - & *fac_shield(i)*fac_shield(j) - ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) - & *fac_shield(i)*fac_shield(j) + endif + ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) + & *fac_shield(i)*fac_shield(j)*sss1 + ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) + & *fac_shield(i)*fac_shield(j)*sss1 C Diagnostics. Comment out or remove after debugging! c ees0p(num_conti,i)=0.5D0*fac3*ees0pij @@ -2343,24 +2550,24 @@ C End diagnostics. c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij, c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont C Angular derivatives of the contact function - ees0pij1=fac3/ees0pij - ees0mij1=fac3/ees0mij - fac3p=-3.0D0*fac3*rrmij - ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij) - ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij) -c ees0mij1=0.0D0 - ecosa1= ees0pij1*( 1.0D0+0.5D0*wij) - ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1) - ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1) - ecosa2= ees0mij1*(-1.0D0+0.5D0*wij) - ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2) - ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2) - ecosap=ecosa1+ecosa2 - ecosbp=ecosb1+ecosb2 - ecosgp=ecosg1+ecosg2 - ecosam=ecosa1-ecosa2 - ecosbm=ecosb1-ecosb2 - ecosgm=ecosg1-ecosg2 + ees0pij1=fac3/ees0pij + ees0mij1=fac3/ees0mij + fac3p=-3.0D0*fac3*rrmij + ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij) + ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij) +c ees0mij1=0.0D0 + ecosa1= ees0pij1*( 1.0D0+0.5D0*wij) + ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1) + ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1) + ecosa2= ees0mij1*(-1.0D0+0.5D0*wij) + ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2) + ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2) + ecosap=ecosa1+ecosa2 + ecosbp=ecosb1+ecosb2 + ecosgp=ecosg1+ecosg2 + ecosam=ecosa1-ecosa2 + ecosbm=ecosb1-ecosb2 + ecosgm=ecosg1-ecosg2 C Diagnostics c ecosap=ecosa1 c ecosbp=ecosb1 @@ -2369,84 +2576,91 @@ c ecosam=0.0D0 c ecosbm=0.0D0 c ecosgm=0.0D0 C End diagnostics - facont_hb(num_conti,i)=fcont - fprimcont=fprimcont/rij + facont_hb(num_conti,i)=fcont + fprimcont=fprimcont/rij cd facont_hb(num_conti,i)=1.0D0 C Following line is for diagnostics. cd fprimcont=0.0D0 - do k=1,3 - dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb) - dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg) - enddo - do k=1,3 - gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k) - gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k) - enddo - gggp(1)=gggp(1)+ees0pijp*xj - gggp(2)=gggp(2)+ees0pijp*yj - gggp(3)=gggp(3)+ees0pijp*zj - gggm(1)=gggm(1)+ees0mijp*xj - gggm(2)=gggm(2)+ees0mijp*yj - gggm(3)=gggm(3)+ees0mijp*zj + do k=1,3 + dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb) + dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg) + enddo + do k=1,3 + gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k) + gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k) + enddo + gggp(1)=gggp(1)+ees0pijp*xj + & +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1 + gggp(2)=gggp(2)+ees0pijp*yj + & +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1 + gggp(3)=gggp(3)+ees0pijp*zj + & +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1 + gggm(1)=gggm(1)+ees0mijp*xj + & +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1 + gggm(2)=gggm(2)+ees0mijp*yj + & +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1 + gggm(3)=gggm(3)+ees0mijp*zj + & +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1 C Derivatives due to the contact function - gacont_hbr(1,num_conti,i)=fprimcont*xj - gacont_hbr(2,num_conti,i)=fprimcont*yj - gacont_hbr(3,num_conti,i)=fprimcont*zj - do k=1,3 + gacont_hbr(1,num_conti,i)=fprimcont*xj + gacont_hbr(2,num_conti,i)=fprimcont*yj + gacont_hbr(3,num_conti,i)=fprimcont*zj + do k=1,3 c c 10/24/08 cgrad and ! comments indicate the parts of the code removed c following the change of gradient-summation algorithm. c cgrad ghalfp=0.5D0*gggp(k) cgrad ghalfm=0.5D0*gggm(k) - gacontp_hb1(k,num_conti,i)=!ghalfp - & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) - & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) - & *fac_shield(i)*fac_shield(j) + gacontp_hb1(k,num_conti,i)=!ghalfp + & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) + & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & *sss1*fac_shield(i)*fac_shield(j) - gacontp_hb2(k,num_conti,i)=!ghalfp - & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) - & *fac_shield(i)*fac_shield(j) + gacontp_hb2(k,num_conti,i)=!ghalfp + & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) + & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & *sss1*fac_shield(i)*fac_shield(j) - gacontp_hb3(k,num_conti,i)=gggp(k) - & *fac_shield(i)*fac_shield(j) + gacontp_hb3(k,num_conti,i)=gggp(k) + & *sss1*fac_shield(i)*fac_shield(j) - gacontm_hb1(k,num_conti,i)=!ghalfm - & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) - & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) - & *fac_shield(i)*fac_shield(j) + gacontm_hb1(k,num_conti,i)=!ghalfm + & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) + & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & *sss1*fac_shield(i)*fac_shield(j) - gacontm_hb2(k,num_conti,i)=!ghalfm - & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) - & *fac_shield(i)*fac_shield(j) + gacontm_hb2(k,num_conti,i)=!ghalfm + & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) + & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & *sss1*fac_shield(i)*fac_shield(j) - gacontm_hb3(k,num_conti,i)=gggm(k) - & *fac_shield(i)*fac_shield(j) + gacontm_hb3(k,num_conti,i)=gggm(k) + & *sss1*fac_shield(i)*fac_shield(j) - enddo - ENDIF ! wcorr - endif ! num_conti.le.maxconts - endif ! fcont.gt.0 - endif ! j.gt.i+1 - if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then - do k=1,4 - do l=1,3 - ghalf=0.5d0*agg(l,k) - aggi(l,k)=aggi(l,k)+ghalf - aggi1(l,k)=aggi1(l,k)+agg(l,k) - aggj(l,k)=aggj(l,k)+ghalf enddo + ENDIF ! wcorr + endif ! num_conti.le.maxconts + endif ! fcont.gt.0 + endif ! j.gt.i+1 +#endif + if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then + do k=1,4 + do l=1,3 + ghalf=0.5d0*agg(l,k) + aggi(l,k)=aggi(l,k)+ghalf + aggi1(l,k)=aggi1(l,k)+agg(l,k) + aggj(l,k)=aggj(l,k)+ghalf + enddo + enddo + if (j.eq.nres-1 .and. i.lt.j-2) then + do k=1,4 + do l=1,3 + aggj1(l,k)=aggj1(l,k)+agg(l,k) enddo - if (j.eq.nres-1 .and. i.lt.j-2) then - do k=1,4 - do l=1,3 - aggj1(l,k)=aggj1(l,k)+agg(l,k) - enddo - enddo - endif - endif + enddo + endif + endif c t_eelecij=t_eelecij+MPI_Wtime()-time00 return end @@ -2455,7 +2669,7 @@ C----------------------------------------------------------------------- C C Compute Evdwpp C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.IOUNITS' @@ -2465,11 +2679,12 @@ C include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' - include 'COMMON.CONTACTS' +c include 'COMMON.CONTACTS' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' - dimension ggg(3) + include "COMMON.SPLITELE" + double precision ggg(3) integer xshift,yshift,zshift c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT @@ -2478,6 +2693,14 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions double precision scal_el /0.5d0/ #endif c write (iout,*) "evdwpp_short" + integer i,j,k,iteli,itelj,num_conti,ind,isubchap + double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb + double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1, + & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, + & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw + double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, + & dist_temp, dist_init,sss_grad + double precision sscale,sscagrad evdw1=0.0D0 C print *,"WCHODZE" c write (iout,*) "iatel_s_vdw",iatel_s_vdw, @@ -2494,12 +2717,12 @@ c call flush(iout) 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.0d0) xmedi=xmedi+boxxsize - ymedi=mod(ymedi,boxysize) - if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize - zmedi=mod(zmedi,boxzsize) - if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i), c & ' ielend',ielend_vdw(i) @@ -2527,44 +2750,44 @@ c call flush(iout) 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 + 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 rrmij=1.0D0/rij rij=dsqrt(rij) c sss=sscale(rij/rpp(iteli,itelj)) c sssgrad=sscagrad(rij/rpp(iteli,itelj)) - sss=sscale(rij) - sssgrad=sscagrad(rij) + sss=sscale(rij/rpp(iteli,itelj),r_cut_respa) + sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa) if (sss.gt.0.0d0) then rmij=1.0D0/rij r3ij=rrmij*rmij @@ -2584,9 +2807,9 @@ 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) + 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) C ggg(1)=facvdw*xj C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj @@ -2617,10 +2840,18 @@ C include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' + include "COMMON.SPLITELE" logical lprint_short common /shortcheck/ lprint_short - dimension ggg(3) + double precision ggg(3) integer xshift,yshift,zshift + integer i,iint,j,k,iteli,itypj,subchap + double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, + & fac,e1,e2,rij + double precision evdw2,evdw2_14,evdwij + double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, + & dist_temp, dist_init + double precision sscale,sscagrad if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb evdw2=0.0D0 evdw2_14=0.0d0 @@ -2635,16 +2866,17 @@ c & ' iatscp_e=',iatscp_e 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)) - 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 + 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) + 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 @@ -2662,14 +2894,14 @@ c corrected by AL zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize c end correction - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize @@ -2681,23 +2913,27 @@ c end correction 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 + 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) - sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) - sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) + sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int) + if (sss1.eq.0) cycle + sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa) + sssgrad= + & sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa) + sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int) if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij), & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss if (sss.lt.1.0d0) then @@ -2707,18 +2943,19 @@ c end correction if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 - evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss) + evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss1 endif evdwij=e1+e2 - evdw2=evdw2+evdwij*(1.0d0-sss) + evdw2=evdw2+evdwij*(1.0d0-sss)*sss1 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') & 'evdw2',i,j,sss,evdwij C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C - fac=-(evdwij+e1)*rrij*(1.0d0-sss) - fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli) + fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1 + fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli) + & +sssgrad1) ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -2762,7 +2999,7 @@ C This subroutine calculates the excluded-volume interaction energy between C peptide-group centers and side chains and its gradient in virtual-bond and C side-chain vectors. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -2773,29 +3010,37 @@ C include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' + include "COMMON.SPLITELE" integer xshift,yshift,zshift logical lprint_short common /shortcheck/ lprint_short - dimension ggg(3) + integer i,iint,j,k,iteli,itypj,subchap + double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, + & fac,e1,e2,rij + double precision evdw2,evdw2_14,evdwij + double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, + & dist_temp, dist_init + double precision ggg(3) + double precision sscale,sscagrad 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,rlamb + if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb do i=iatscp_s,iatscp_e 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)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) - 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 + 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 c if (lprint_short) c & write (iout,*) "i",i," itype",itype(i),itype(i+1), @@ -2803,7 +3048,7 @@ c & " nscp_gr",nscp_gr(i) do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) c if (lprint_short) c & write (iout,*) "j",j," itypj",itypj if (itypj.eq.ntyp1) cycle @@ -2823,18 +3068,18 @@ c corrected by AL zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize c end correction - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 c if (lprint_short) then c write (iout,*) i,j,xi,yi,zi,xj,yj,zj c write (iout,*) "dist_init",dsqrt(dist_init) c endif - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize @@ -2846,24 +3091,25 @@ c endif zj_temp=zj subchap=1 endif - enddo - enddo - enddo + enddo + enddo + enddo c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp) - 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 + 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=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) - sss=sscale(1.0d0/(dsqrt(rrij))) - sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) + sss=sscale(1.0d0/(dsqrt(rrij*rscp(itypj,iteli))),r_cut_respa) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij*rscp(itypj,iteli))), + & r_cut_respa) if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij), & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),